Molecular Dynamics for Peptide Conformational Sampling

A molecular dynamics trajectory is not a measurement. It is the output of integrating equations of motion under an empirical potential, and the conformational population it produces is that potential's prediction about your peptide, sampled for a finite time. Two independent things can be wrong with it: the potential can be inaccurate, and the sampling can be incomplete. Convergence analysis addresses only the second: a best-practices reference on uncertainty in molecular simulation covers trajectory data from conventional molecular dynamics and single-chain Monte Carlo, offers guidance for some enhanced sampling approaches, and states that it does not discuss systematic errors arising, for example, from inaccuracy in the chosen model or force field ^[1].

This is a different question from the one docking asks. Docking predicts a bound geometry and ranks it; conformational sampling asks what fraction of time a peptide occupies each region of its conformational space, and what free-energy cost separates those regions.

The force field is part of the result

In a benchmark of 21 systems spanning folded and disordered proteins, simulated with six force fields and compared against more than 9,000 experimental data points, none of the tested force fields simultaneously provided accurate descriptions of folded proteins, of the dimensions of disordered proteins, and of the secondary structure propensities of disordered proteins ^[2]. That is a finding about those six, not a permanent ceiling: guided by results on a subset of the benchmark, the same authors modified one force field and report that the result, a99SB-disp, achieved excellent agreement with experiment for disordered proteins while maintaining state-of-the-art accuracy for folded ones ^[2].

The magnitude of the force-field effect can be larger than expected. Eight all-atom force fields were compared on intrinsically disordered protein ensembles against primary small-angle X-ray scattering and NMR data; the ensembles differed markedly in chain dimensions, hydrogen bonding and secondary structure content, and changing the force field had a stronger effect on secondary structure content than changing the entire peptide sequence ^[3]. A high population of left-handed α-helix appeared in the CHARMM 36 ensemble, inconsistent with measured scalar couplings, while the CHARMM 22* ensemble had the lowest error in chemical shifts and J-couplings and agreed well with the scattering data ^[3].

Part of the error in disordered-state dimensions lives in the water model. Writing in 2015, the developers of TIP4P-D observed that MD simulations using then-current physics-based models tend to produce disordered-state ensembles that are structurally too compact relative to experiment; they found that the water models typically used in MD significantly underestimate London dispersion interactions, and speculated that this may be a possible reason for those results ^[4]. Simulations of solvated proteins with the corrected water model typically produced disordered states that were substantially more expanded and in better agreement with experiment ^[4]. A radius of gyration is therefore a joint property of protein force field and water model, and reporting one without the other omits what another group would need in order to reproduce it.

Helical peptides fail differently. Backbone energy corrections were determined for two force fields so as to reproduce the fraction of helix measured in short peptides at 300 K, and the optimised force fields gave results in excellent agreement with NMR for folded proteins and for short peptides not used in the optimisation ^[5]. Agreement at one temperature did not repair the thermodynamics: the dependence of helix content on temperature was too weak, which the authors describe as a problem shared with other force fields, and a fit of Lifson-Roig helix-coil theory showed the helix extension parameter w agreeing well with experiment while its entropic and enthalpic components were both only about half the respective experimental estimates ^[5]. A temperature-dependent helical fraction can therefore still be wrong even when the 300 K helicity looks correct.

Explicit versus implicit solvent

The generalized Born model is one of the fastest implicit solvent models and is widely adopted for MD ^[6]. A comparison of explicit-solvent particle mesh Ewald with TIP3P against a popular generalized Born model, both in AMBER, found conformational sampling speedups that were highly system- and problem-dependent: at matched temperature, approximately onefold for small changes taken as dihedral angle flips in a protein, between about 1- and 100-fold for large changes taken as nucleosome tail collapse and DNA unwrapping, and approximately sevenfold for a mixed case, the folding of a miniprotein ^[7]. Only the mixed case is peptide-scale — CLN025, which the authors describe as the smallest fast-folding miniprotein at 10 residues; the small and large cases are a roughly 300-residue protein and a nucleosome.

Two further results from that work matter more than the speedups. For all systems studied, the sampling speedup was mainly due to reduction in solvent viscosity rather than possible differences in free-energy landscapes between the solvent models; and those landscapes may differ substantially between the models ^[7]. If the acceleration is largely a viscosity effect, a relaxation time read off an implicit-solvent trajectory is not the same quantity as one read off an explicit-solvent trajectory, and comparing either to an experimental rate takes an argument rather than an assumption.

Secondary-structure bias in implicit solvent is a parameterisation property, not a fixed defect of the approach. The authors of GB-Neck2 note that the speed of generalized Born models comes with tradeoffs and that many reports have pointed out weaknesses with them; they refit the empirical parameters and report that the resulting model's secondary structure preferences are in much better agreement with those from explicit-solvent MD, alongside near-quantitative reproduction of experimental structure and thermal stability profiles for several model peptides with varying secondary structure motifs ^[6]. The operative question is which variant you are running, and whether it was validated on peptides of your class.

The timescale gap

Atomic-level simulations over periods ranging between 100 μs and 1 ms, conducted with a single physics-based energy function, showed 12 structurally diverse proteins folding spontaneously and repeatedly to their experimentally determined native structures ^[8]. That is the benchmark against which to judge whether your process is accessible to unbiased simulation at all.

Some peptide processes are. For the short cyclic peptide Ac-(cyclo-2,6)-R[KAAAD]-NH₂, microsecond-timescale MD analysed by Cartesian-coordinate principal component analysis, with a Markov model built from that analysis, gave a free energy of unfolding of −0.57 kJ/mol and a relaxation time of about 0.435 μs at 298.15 K, against experimentally reported values of −0.22 kJ/mol and 0.42 μs ^[9]. A single-turn helix in a constrained cyclic peptide is close to the best case: a fast process in a small, conformationally restricted system. The question worth asking is not how long the trajectory ran, but how many times it crossed the barrier you care about.

Enhanced sampling, and what each surface establishes

Replica exchange

For systems whose long-time dynamics is dominated by slow interconversion between two metastable states, the relative efficiency of replica exchange molecular dynamics and plain MD at comparable computational cost is given by the ratio of the number of transitions between the two states averaged over all replicas at the different temperatures to the number of transitions at the single temperature of the MD run, provided exchange is frequent compared with the transition times ^[10]. High efficiency is achieved by including replica temperatures at which transitions are more frequent than at the temperature of interest ^[10]. A ladder whose top temperature never crosses the barrier buys parallel replicas and no extra sampling — though that derivation covers two-state systems with frequent exchange, not arbitrarily rugged landscapes.

Solute tempering changes the arithmetic: in a direct comparison the number of replicas could be reduced four to five times relative to conventional replica exchange without affecting the temperature range or compromising an efficient random walk of replicas over temperatures ^[11]. The convergence result from that study is what to carry into peptide work. Solute tempering reproduced conventional replica exchange data extremely well for two Aβ_10-40 peptides binding a zwitterionic dimyristoylphosphatidylcholine bilayer, while agreement for the short Aβ_25-35 fragment in water was less perfect and improved with additional solute-tempering simulations; the authors call it surprising that convergence was much better for the ordered peptides on the bilayer than for the small unstructured peptide in water ^[11]. A small disordered peptide is not an easy sampling problem merely because it is small.

Metadynamics

Metadynamics explores a multidimensional free energy surface in the space defined by a few collective coordinates using a history-dependent potential term that, in time, fills the minima in the surface, allowing efficient exploration and accurate determination of the surface as a function of those coordinates ^[12]. Well-tempered metadynamics determines that dependence using an adaptive bias, giving a formalism with metadynamics and canonical sampling as limiting cases, in which convergence and errors can be rigorously and easily controlled ^[13].

Collective-variable choice is one place free-energy surfaces go wrong, and the failure is silent. Umbrella sampling, adaptive biasing force and metadynamics all bias predefined, system-dependent collective variables, and a review of the alternatives characterises precisely this class as subject to the hidden energy barrier problem ^[14]. A surface that looks smooth and well sampled along the variable you biased can conceal a barrier in a degree of freedom you did not. That review groups replica exchange and parallel tempering among methods that do not require predefined system-dependent collective variables and as such do not suffer from that problem ^[14] — a reason to run both classes rather than trust either alone.

Umbrella sampling

Umbrella sampling biases individual windows along a set of order parameters and recombines their biased probability densities. The weighted histogram analysis method is the standard protocol for postprocessing that information into a potential of mean force with respect to those order parameters, determining the unbiased density of states by satisfying a self-consistent condition through an iterative procedure ^[15]. It works very effectively when the number of order parameters is small, and its computational cost grows rapidly in higher dimension — a limitation that motivated a multivariate linear regression alternative shown to give landscapes comparable in accuracy at a small fraction of the cost ^[15]. A potential of mean force is a projection: the free energy along the order parameters you chose, silent about barriers orthogonal to them.

What a converged surface does not establish

Equilibrium simulations of a fast-folding villin headpiece variant with four different force fields each produced a large number of transitions between unfolded and folded states, and in all four cases both the folding rate and the native structure agreed well with experiment. The folding mechanism and the properties of the unfolded state nonetheless depended substantially on the choice of force field, and the authors conclude that although it is important to match a single, experimentally determined structure and folding rate, this does not ensure that a given simulation will provide a unique and correct description of the full free-energy surface and the mechanism of folding ^[16].

Convergence, and why one long trajectory is not evidence

A single trajectory is one realisation of a stochastic process, and a plateau in whatever observable you happened to plot is weak evidence of equilibrium, because you chose the projection. The expectation set out in the best-practices literature is that modellers analyse and communicate statistical uncertainties so that consumers of simulated data understand its significance and limitations ^[1]. The apparatus is independent runs from independently prepared starting structures, decorrelation-time estimates for the observables reported, and error bars from between-run variance rather than within-run point counts.

The eight-force-field comparison shows that standard applied properly: sampling uncertainty was bounded rather than assumed away, with 0.964 ms of aggregate simulation used to show that residual sampling uncertainty was much smaller than the between-force-field differences being reported ^[3]. Where a kinetic claim is made, a discretised state model makes it falsifiable, as in the cyclic peptide study whose Markov model produced the relaxation time then compared with experiment ^[9]. None of this touches the potential: uncertainty analysis explicitly excludes systematic errors arising from inaccuracy in the chosen model or force field ^[1].

Validating against experimental observables

Chemical shifts and scalar couplings

Force fields have been evaluated against a suite of 524 chemical shift and J coupling measurements on dipeptides, tripeptides, tetra-alanine and ubiquitin; two of the force fields examined achieved high accuracy, and for those two the calculation error was comparable to the uncertainty in the experimental comparison — which the authors say suggests that extracting further force field improvements from NMR data may require increased accuracy in J coupling and chemical shift prediction ^[17]. That defines the floor you work against. SPARTA+ reports standard deviations between predicted and experimental shifts of 2.45, 1.09, 0.94, 1.14, 0.25 and 0.49 ppm for δ¹⁵N, δ¹³C′, δ¹³C^α, δ¹³C^β, δ¹H^α and δ¹H^N respectively, on a set of eleven validation proteins, and its authors describe it as a modest but consistent improvement of 2–10% over the best programs then available that appears to approach the limit at which empirical approaches can predict chemical shifts ^[18]. Those figures come from validation proteins, not peptides, so treat them as an order-of-magnitude floor rather than a peptide-specific tolerance. Reproducing a shift to better than the predictor's own scatter is not evidence the ensemble is right.

Order parameters

Apparent systematic discrepancies between order parameters extracted from simulations and from experiments are common, particularly for elements of noncanonical secondary structure ^[19]. A 1.2 μs explicit-solvent simulation of ubiquitin compared against backbone order parameters previously derived from NMR relaxation revealed fluctuations in three loop regions occurring on timescales comparable to or longer than the overall rotational diffusion of ubiquitin, whose effects would not be apparent in the experimentally derived order parameters; a coupled analysis of internal and overall motion yielded simulated order parameters substantially closer to the experimental values than a conventional analysis of internal motion alone ^[19]. The procedural point transfers even though that system is a protein: before attributing an S² discrepancy to the force field, establish whether the motion is comparable to or slower than overall rotational diffusion, because motion on that timescale would not be apparent in the experimental order parameter.

Circular dichroism

A CD spectrum can be computed from a structural model or ensemble using the average secondary structure composition together with a precalculated set of basis spectra, an approach designed for model validation — to estimate the error of a structural model against a measured CD spectrum ^[20]. For over 80% of the globular reference proteins in that work the basis sets accurately predicted the experimental spectrum solely from secondary structure composition; for the remaining 20%, correcting for intensity normalisation considerably improved prediction power ^[20]. For short peptides, and for an example complex of intrinsically disordered proteins, predictions strongly benefited from accounting for side-chain contributions and structural flexibility ^[20]. A peptide CD spectrum back-calculated from secondary-structure content alone is the weakest form of this comparison, and a mismatch under that treatment does not by itself convict the ensemble.

Why agreement is weak evidence

Structural ensembles often offer more useful representations of proteins than individual conformations, and the principles of protein structure determination have to be extended to handle conformationally heterogeneous states and experimental measurements that are averaged over such states and affected by a variety of errors ^[21]. Because the comparison sets one experimental number against an average over a whole ensemble, many different ensembles reproduce it — and objective criteria for the comparison, validation, visualisation and dissemination of ensembles are described in that review as a future goal rather than settled practice ^[21]. Agreement with an ensemble-averaged observable is consistent with your ensemble; it does not single it out.

What to report

• Force field version, water model and ion parameters, and whether that pairing was validated on systems resembling yours.
• Aggregate simulation time, number of independent runs, and how their starting structures were generated.
• The observable used to judge convergence, its decorrelation time, and the between-run spread rather than a single-run standard error.
• For metadynamics: collective variables, bias deposition parameters, and evidence that an orthogonal slow degree of freedom was checked. For replica exchange: temperature range, exchange acceptance rate, round trips through the ladder, and transitions across the barrier of interest. For umbrella sampling: order parameters, window spacing, restraint force constants and window overlap.
• For every back-calculated observable: the predictor used and its published error, stated next to the agreement claimed.