Molecular Docking and In Silico Prediction for Peptides: What a Score Is Not

Molecular docking asks a structural question — given a receptor structure and a ligand, what is the most plausible bound geometry? — and returns a number, the score a scoring function assigns to that geometry. The two are often conflated, and that is where most misreadings of a docking result begin.

This guide covers why peptide flexibility violates the assumptions small-molecule docking engines were built on, what published peptide docking benchmarks actually report, how pose quality is measured, what a docking score is and is not, and what it takes to move a prediction from hypothesis to supported claim.

Why peptides break the small-molecule assumption

Small-molecule docking engines treat the receptor as rigid or nearly so, enumerate ligand conformations by rotating a modest number of torsions, and score each placement. That works when the ligand has a handful of rotatable bonds and a preferred low-energy shape it mostly keeps. Peptides do not fit that description. The authors of one benchmarking study describe docking studies on protein-peptide interactions as challenging and time-consuming because peptides are generally more flexible than proteins and tend to adopt numerous conformations, and note that a series of docking methods had not been rigorously validated for protein-peptide complexes in the literature, unlike the protein-protein and protein-ligand cases ^[1]. Another assessment describes even short peptides as difficult targets owing to their high torsional flexibility ^[2]. And the conformation matters twice over: peptides often adopt different conformations in their free and bound forms, so a single starting structure encodes an assumption about which one you are docking ^[3].

A coarse combinatorial illustration makes the scale visible. Suppose, purely as a supposition, that each residue contributes two independently variable backbone dihedrals and each dihedral is discretized into three states. A 10-residue peptide then has 2 x 10 = 20 dihedrals and 3^20 = 3,486,784,401 backbone states, before any side chain. A small molecule with 6 rotatable bonds, discretized the same way, has 3^6 = 729. The ratio is 3^14 = 4,782,969 — roughly a five-million-fold difference in the space a search must cover, from backbone alone. This arithmetic is illustrative and uncited. It overstates the problem, because backbone dihedrals are correlated and large regions of the space are sterically forbidden, and understates it, because real dihedrals are continuous and side chains add their own.

Two design responses follow from this. One is to fold and dock at once: AutoDock CrankPep folds a peptide in the potential field created by the protein rather than docking a pre-built conformer, an approach its authors report on the LEADS-PEP and peptiDB datasets of peptides up to 15 residues and extend to peptides of 16 to 20 residues ^[4]. The other is to dock an ensemble: a HADDOCK protocol combines the most frequently found peptide conformations in complexes with representative conformations drawn from molecular dynamics of the free peptide, and its authors demonstrate it on six cases of increasing difficulty taken from a protein-peptide benchmark set; in the majority of cases where MD conformations were added to the starting ensemble, they report an improvement in the quality of the resulting models ^[3]. Note that this second protocol is information-driven: knowledge of the binding site on the receptor drove the docking in each case ^[3], a constraint that should always be stated.

What peptide docking benchmarks actually report

Benchmarks are the only honest source of expectation here, and the published ones are sobering. A benchmarking study of six docking methods on 133 protein-peptide complexes with peptides 9 to 15 residues long evaluated results with the CAPRI parameters FNAT, I-RMSD and L-RMSD. In blind docking, the best-performing method by top pose was FRODOCK, with an average L-RMSD of 12.46 Å. In re-docking, ZDOCK gave an average L-RMSD of 8.60 Å for the top pose ^[1]. An average ligand RMSD in that range is not a near-miss on the binding mode; it is a different answer.

A separate assessment compiled 47 complexes with peptides of up to five residues from 11 studies and benchmarked AutoDock Vina against them, using a deliberately stringent measure: the quality of the top-scoring pose rather than the best pose anywhere in the output ^[2]. Those authors are explicit that the varying strategies and constraints of the source studies preclude direct quantitative comparison between them, and they note that guidelines on how much sampling is required for a reproducible result were, at the time, lacking ^[2]. If the sampling required for reproducibility is unestablished, two runs of the same protocol on the same input can disagree, and a single run reports the search as much as it reports the system.

Structure prediction networks have moved these numbers. On a set of 112 peptide-protein complexes, AlphaFold-Multimer produced models of acceptable or better quality (DockQ at or above 0.23) for 66 complexes, 25 of them high quality (DockQ at or above 0.8) — against 23 or 47 acceptable and only 4 or 8 high quality for energy-based docking and template-based approaches respectively ^[6]. That is 58.9 percent acceptable or better and 22.3 percent high quality — a large improvement on what came before, and still roughly two in five complexes with no acceptable model.

Sampling failure and ranking failure are different failures

The most useful diagnostic is the gap between the top-ranked pose and the best pose the run produced. In the 133-complex study, FRODOCK's blind docking average L-RMSD was 12.46 Å for the top pose but 3.72 Å for the best pose — a gap of 8.74 Å. For ZDOCK re-docking the figures were 8.60 Å and 2.88 Å, a gap of 5.72 Å ^[1]. Those authors conclude directly that the ranking of generated poses needs improvement and that the present ranking scheme is not satisfactory ^[1].

The same split appears in the network results. Perturbing the AlphaFold network weights to force broader sampling raised the number of acceptable models from 66 to 75 of 112 and the median first-ranked DockQ from 0.47 to 0.55, while the best possible DockQ rose from 0.58 to 0.72; the authors note that selecting the best possible model remains a challenge ^[6]. A first-ranked median of 0.55 against a best possible DockQ of 0.72 is a selection problem, not a search problem.

This changes what a negative result means. If the correct pose is in your output but ranked eleventh, the search worked and the scoring failed. If it is absent, no amount of rescoring helps. A protocol reporting only the top pose cannot tell these apart.

How pose quality is measured, and why one number hides three

The CAPRI convention assesses a docking model with three related but independent quantities: FNAT, the fraction of native contacts recovered; L-RMSD, the ligand RMSD after superposing the receptor; and I-RMSD, the interface RMSD. Models are binned as incorrect, acceptable, medium or high quality by applying cutoffs to these three ^[5]. The quantities can disagree. A model with a relatively poor L-RMSD above 10 Å can still qualify as acceptable given a decent FNAT above 0.50 and an I-RMSD below 3.0 Å ^[5]. DockQ combines the three into a single continuous score on the interval 0 to 1; its authors report it reproduces the CAPRI classification with an average positive predictive value of 94 percent at 90 percent recall, and they describe the CAPRI cutoffs themselves as ad hoc ^[5].

The practical consequence: an unqualified RMSD is ambiguous until you say which RMSD, after which superposition, over which atoms. Two reports of 2.5 Å may not be the same quantity.

What a docking score is not

The clearest evidence comes from the small-molecule side, where benchmarks are largest. CASF-2016 evaluates scoring functions on 285 protein-ligand complexes with high-quality crystal structures and reliable binding constants, deliberately decoupling scoring from docking, and grades each function on four separate metrics it calls scoring power, ranking power, docking power and screening power. Testing 25 scoring functions, its authors report that current performance is more promising in terms of docking power than scoring, ranking or screening power ^[7]. In the field's standard usage — stated here without citation — docking power is picking a near-native pose out of decoys, scoring power the correlation between predicted and measured affinity, ranking power the correct ordering of ligands for one target, and screening power the enrichment of true binders from a library.

The task these functions do best is the geometric one. It is not the task people quote them for. A docking score is not an affinity, not a free energy, and not comparable across different receptors; it is a pose-selection statistic reported in energy-like units.

The units invite a conversion that is worth doing once, because it sets the precision bar. At 298.15 K, RT is 8.314 x 298.15 = 2,479 J/mol, or 0.592 kcal/mol. A tenfold change in a dissociation constant corresponds to RT ln(10) = 0.592 x 2.303 = 1.36 kcal/mol. So distinguishing a 100 nM binder from a 1 µM binder requires a method whose error is meaningfully below 1.36 kcal/mol. Docking scores are not characterized to that standard. This conversion is standard thermodynamics, and the comparison it supports — between the precision an affinity claim requires and what a pose-selection statistic offers — is drawn here rather than made by any of the cited sources; both are stated without citation.

What scoring functions learn when the benchmark is biased

Machine-learned scoring functions raise a distinct problem. A re-evaluation of three widely used structure-based virtual screening benchmarks concluded that not every benchmark dataset is suitable for machine learning, and demonstrated on two examples from the literature that bias is learned implicitly and unnoticed from standard benchmarks ^[8]. One evaluation of Gnina, a docking program using deep convolutional networks for scoring, found its default scoring outperformed the empirical AutoDock Vina scoring function on 89 of 117 DUD-E and LIT-PCBA targets (76.1 percent), with a median 1 percent early enrichment factor more than twice Vina's, and stated that bias lingers in these sets even when not used directly to train models, obfuscating the extent to which performance comes from interpreting molecular interactions rather than fitting non-informative property distributions ^[9]. A separate group reported that models trained on their modeled-complex dataset could reduce the bias caused by buried solvent-accessible surface area that they had previously found in models trained on PDBbind ^[10].

These are small-molecule benchmarks, not peptide ones, which is part of the point: peptide docking protocols routinely borrow scoring functions fitted on small-molecule complexes, inheriting whatever those functions actually learned.

The entropy term the score mostly does not have

A score computed on a single rigid bound pose contains no representation of the conformational freedom the free peptide gave up in order to adopt it — that is a consequence of the calculation's inputs, not a finding, and it is stated here without citation. The magnitude of what is missing has been estimated. Studying the N-terminal SH3 domain of CrkII with a proline-rich motif derived from C3G, and combining NMR model-free analysis with delta-2D, DynaMine and a structure-based entropy calculation, one group concluded that the contribution of backbone conformational entropy change is significant in interactions mediated by intrinsically disordered proteins and regions, which undergo a disorder-to-order transition on binding that carries a large conformational entropy penalty ^[11]. Those authors reached this from a combination of methods on one domain-motif pair; it is not a general constant, and a score computed on a single rigid pose carries no term for it.

Rescoring and free-energy methods do better, and are still not free

Physics-based free-energy methods narrow the gap. A free-energy perturbation study whose retrospective results encompass 200 ligands framed the target explicitly: the accuracy needed to guide lead optimization is about fivefold discrimination in binding affinity, which its authors describe as long challenging to achieve reliably across diverse ligands and targets ^[12]. Fivefold in affinity is 1.36 x log10(5) = 0.95 kcal/mol — under one kilocalorie per mole.

For a protein-protein case, an FEP protocol applied to antibody-antigen complexes was tested against three surface plasmon resonance datasets in which antibody residues in the interface were systematically mutated to alanine. Across 55 cases the RMS error versus experiment was 0.68 kcal/mol, which the authors call near experimental accuracy; the correlation coefficient for the combined dataset including residues with glycan contacts was R2 = 0.49, which they judge sufficient to guide residue selection, with the stated caveat that this assumes the accuracy can be realized prospectively ^[13]. An RMS error of 0.68 kcal/mol corresponds to 10^(0.68/1.36) = 3.2-fold in affinity. That is a computationally expensive, carefully executed method on a well-characterized system, and the residual uncertainty is roughly three-fold, not three percent.

How a prediction gets validated

Nothing in the sections above validates a docking result. Validation requires a measurement the prediction did not see, and the following practices are convention, not findings from any single primary study.

• Re-docking a ligand into the crystal structure it came from is the weakest available check. The receptor is already shaped to the answer.
• Report sampling explicitly: number of independent runs, random seeds, exhaustiveness or equivalent, and whether repeat runs converged. Given that guidelines on the sampling required for reproducibility were reported as lacking when peptide docking was assessed against that question ^[2], a single run is an anecdote.
• Report the pose ranked first and the best pose generated, with the quality metric named (FNAT, L-RMSD, I-RMSD or DockQ). The gap between them is the diagnostic ^{[1, 5]}.
• Report every constraint imposed. An information-driven run that used a known binding site is answering a narrower question than a blind run ^[3].
• Test the model with an orthogonal experiment that can falsify it. Interface mutagenesis with binding measured by surface plasmon resonance is the pattern used in the FEP study above ^[13]: the prediction makes a per-residue claim, and the mutant either shifts the measured affinity or it does not.

What is not known

As of that 2015 assessment, guidelines on how much sampling a flexible peptide docking run needs before its result is reproducible were lacking, and its authors presented their work as the first study to address the question in detail ^[2]. How to rank the poses a search produces remains unsolved — the explicit conclusion of the 133-complex benchmark ^[1], and visible in the residual gap between first-ranked and best-possible models from structure prediction networks ^[6]. Whether the gains of machine-learned scoring functions reflect learned interaction physics or benchmark artifacts is stated as unresolved by the authors of the evaluation that reported those gains ^[9]. And whether free-energy accuracy demonstrated retrospectively transfers to prospective use is flagged as an assumption by the authors who reported it ^[13].

A docking model is a structural hypothesis with a provenance. Reported with its search parameters, constraints, quality metric and rank gap, it is a useful one. Reported as a single score in kilocalories per mole, it is a number whose meaning has been discarded.