Protein Descriptors for Machine Learning in Drug Discovery

When you’re modeling protein-ligand binding or reactivity, one of the first decisions is what to actually compute as a protein feature. The options range from a single vector derived from the amino acid sequence to a full ensemble of structural and interaction features extracted from molecular dynamics. Picking the wrong level of detail either throws away information the model needs, or buries it in noise the model can’t use.

The descriptors below split into three fundamentally different categories depending on what inputs they require and how they respond to target movement:

• Invariant (sequence-only): Computed once per protein, identical across every conformation.
• State-dependent (structure/MD): Change from one snapshot to the next, and only exist once you have a structure.
• Complex-dependent (interaction): Require a bound ligand complex, changing with conformation and ligand pose.

The practical recipe that falls out of this split: use sequence descriptors as your base feature layer, then add per-snapshot structural and interaction features on top, aggregated across the MD ensemble (mean, std, percentiles, occupancies) rather than taken from a single frame.

1. Sequence-based descriptors (invariant across MD snapshots)

These only need the amino acid sequence, so they’re cheap and don’t change no matter which conformation you’re looking at.

• AAindex / Z-scales: Physicochemical traits per residue (hydrophobicity, polarity, volume), compressed into a small set of scales.
• BioVec / ProtVec: Embeddings that capture local sequence context and motifs, similar in spirit to word2vec.
• PSSM: Evolutionary conservation pulled from multiple target sequence alignments.
• AA Composition (AAComp), CTD (CTDC/CTDT/CTDD), PAAC: Classical composition and order-aware descriptors, covered in detail below.
• Deep protein language models (PLMs): UniRep, ESM, ProtBERT/ProtT5, and AlphaFold2 embeddings. These are dense numerical vectors encoding biochemical and structural information straight from the raw sequence.

Protein LM embeddings in particular are useful for:

• Cross-species and cross-target generalization.
• Substrate/inhibitor specificity classification.
• Transfer learning for low-data enzymes.

Mutation indicators, a related but separate idea, are binary vectors flagging known catalytic or binding-site mutations. They’re also invariant per variant, so they belong in this same base layer.

2. Structure-based descriptors (change with conformation)

These need at least one 3D structure, and they shift from snapshot to snapshot in an MD trajectory.

• Secondary structure content, especially in loops around the active site.
• Backbone $\phi$/$\psi$ angle distributions.
• Distances describing access channels and gates.
• Pocket geometry: volume and bottleneck radius.

3. Complex-based descriptors (need a protein-ligand structure)

These require a bound complex, not just a protein structure on its own.

• Interaction fingerprints (IFPs): H-bonds, $\pi$-stacking, cation-$\pi$, hydrophobics, and salt bridges (e.g., PLIP, BINANA, ODDT).
• PLEC fingerprints (protein-ligand extended connectivity): Encode local atom environments for protein-ligand contacts within a cutoff. Useful for virtual screening, affinity prediction, and QSAR on docked or MD complexes, but not applicable to sequence-only modeling.

[!NOTE] A related idea worth flagging: Prot2Drug, a transformer decoder that generates ligands directly from protein sequence embeddings. It is validation that sequence-only drug-target interaction (DTI) pipelines can generalize to unseen proteins.

Propy3: a practical toolkit for sequence descriptors

Propy3 computes the classical sequence descriptor families referenced above. Here’s what each one actually encodes:

• Amino Acid Composition (AAComp) - 20 descriptors: Fraction of each amino acid in the sequence. Global composition, no order information.
• CTD descriptors (Composition-Transition-Distribution) - 147 descriptors:
  • CTDC (Composition) - 21: Physicochemical class composition (hydrophobicity, volume, polarity, etc.).
  • CTDT (Transition) - 21: Transition frequencies between classes, adding simple order information.
  • CTDD (Distribution) - 105: Positions of class occurrences along the sequence (1/25/50/75/100%).
• Pseudo amino acid composition (PAAC) - 30 descriptors: Adds correlation factors on top of composition to capture local sequence-order effects.

Which of these to actually keep for binding modeling

Not all five families pull equal weight. Here’s the practical call:

PAAC and the CTD composition/transition families are the strongest sequence-only baseline. CTDD and raw AAComp are optional additions, worth including if you have the budget but not essential.

Layering in the state-dependent features: MD / ensemble descriptors

Sequence descriptors (AAComp/CTD/PAAC) won’t move across MD snapshots, by construction. Everything that actually responds to conformation has to be computed per snapshot and then aggregated over the ensemble. For binding-relevant modeling, that means:

Active-site geometry

• Distance from ligand reaction site to catalytic center (e.g. metal atom or active site residue).
• Local pocket residues planarity; substrate centroid to active site center distance.
• Pocket volume / bottleneck radius (computed with tools like CAVER).

Access and flexibility

• Loop opening distances around the pocket.
• RMSF per residue (loops around active site).
• Secondary structure fractions (DSSP) in those regions.

Interactions (complex-aware)

• IFPs: H-bonds, $\pi$-stacking, cation-$\pi$, hydrophobics, and salt bridges (using PLIP/BINANA/ODDT).
• PLEC fingerprints computed per snapshot.
• Contact counts with key residues; water-mediated contacts.

Surface/electrostatics

• SASA (protein, pocket, substrate, reactive atom).
• Pocket hydrophobicity; APBS potential at substrate.

Dynamics summaries

• % of frames meeting a catalytic cutoff (e.g., reaction site distance $< 4.5\text{ \AA}$).
• Dwell times in "open/closed" channel states.
• MSM state probabilities (optional).

A practical feature blueprint (one row per simulation/complex)

1. Data inputs:
   • Structure/MD: top.pdb, traj.xtc (or a single PDB).
   • Selections: Active center atom (e.g., catalytic metal, or pocket centroid CA), ligand LIG, pocket residues (within 6-8 Å of the pocket).
2. Feature blocks:
   • 3D shape (pocket-level): 3D Zernike moments on the pocket surface, aggregated over MD (mean $\pm$ std).
   • Geometry (explicit metrics): Distance from active center to reacting atom, ligand centroid to pocket reference plane, pocket volume/bottleneck, SASA, aggregated as mean/std/min/max plus % of frames meeting cutoffs.
   • Interactions: PLEC (cutoff ~4.5-5.0 Å) plus PLIP/BINANA contact types, aggregated by occupancy.
3. Output: A single feature table (CSV/parquet) combining continuous descriptors with fingerprint occupancies.
4. Practical cutoffs and tips:
   • Analyze every 10-50 ps (thin the trajectory).
   • Normalize continuous features ($z$-score).
   • Reduce high-dimensional shape features with PCA.

Code: computing basic geometry features with MDAnalysis

Here is a minimal python starting point for computing reference distance and ligand-to-pocket-plane distance, thinned to every 10th frame:

import MDAnalysis as mda
import numpy as np
import pandas as pd

u = mda.Universe("top.pdb", "traj.xtc")

# Select reference atom in the active site (e.g., pocket center CA)
ref_atom = u.select_atoms("protein and resid 102 and name CA")[0]
lig = u.select_atoms("resname LIG")
# Select specific reacting ligand atom
react_atom = lig.select_atoms("name C1")[0]

# Select pocket residues defining a local plane
plane_atoms = u.select_atoms("protein and name CA and (resid 90 or resid 95 or resid 100)")

rows = []
for ts in u.trajectory[::10]:  # thin trajectory
    # 1. Reference distance
    d_ref = np.linalg.norm(ref_atom.position - react_atom.position)
    
    # 2. Ligand-to-pocket-plane distance
    X = plane_atoms.positions - plane_atoms.positions.mean(0)
    _, _, vh = np.linalg.svd(X, full_matrices=False)
    n = vh[-1]  # normal vector of the plane
    p0 = plane_atoms.positions.mean(0)
    
    lc = lig.center_of_mass()
    d_plane = abs(np.dot(lc - p0, n))
    
    rows.append({
        "frame": ts.frame, 
        "d_ref": d_ref, 
        "d_plane": d_plane
    })

pd.DataFrame(rows).to_csv("md_geom.csv", index=False)

Shape descriptors: Zernike moments vs ReOmokage feature points

The choice comes down to what you’re doing downstream: Zernike moments feed cleanly into an ML pipeline, while ReOmokage-style landmarks are better suited to structural alignment tasks where the actual spatial arrangement matters.

Key Takeaways

• Treat sequence descriptors (PAAC, CTDC, CTDT as the core; CTDD, AAComp as optional) as the invariant base layer of any protein feature set.
• Layer structure- and interaction-based features on top only when you have a structure or a bound complex, and compute them per MD snapshot.
• Never take structural or interaction features from a single frame. Aggregate over the ensemble (mean/std/percentiles/occupancies), since a single snapshot is a poor proxy for a dynamic active site.
• PLEC and IFP-style interaction fingerprints are the highest-value additions when you have MD trajectories of a bound complex, since they capture exactly the contacts that a docking pose or a single crystal structure would miss.