Stochastic Path Approach to Compute Atomically Detailed Trajectories: Application to the Folding of C Peptide

Neural Network for Principle of Least Action

The principle of least action is the cornerstone of classical mechanics, theory of relativity, quantum mechanics, and thermodynamics. Here, we describe how a neural network (NN) learns to find the trajectory for a Lennard-Jones (LJ) system that maintains balance in minimizing the Onsager–Machlup (OM) action and maintaining the energy conservation. The phase-space trajectory thus calculated is in excellent agreement with the corresponding results from the “ground-truth” molecular dynamics (MD) simulation. Furthermore, we show that the NN can easily find structural transformation pathways for LJ clusters, for example, the basin-hopping transformation of an LJ₃₈ from an incomplete Mackay icosahedron to a truncated face-centered cubic octahedron. Unlike MD, the NN computes atomic trajectories over the entire temporal domain in one fell swoop, and the NN time step is a factor of 20 larger than the MD time step. The NN approach to OM action is quite general and can be adapted to model morphometrics in a variety of applications.

A survey of multiscale modeling: Foundations, historical milestones, current status, and future prospects

Research problems in the domains of physical, engineering, biological sciences often span multiple time and length scales, owing to the complexity of information transfer underlying mechanisms. Multiscale modeling (MSM) and high-performance computing (HPC) have emerged as indispensable tools for tackling such complex problems. We review the foundations, historical developments, and current paradigms in MSM. A paradigm shift in MSM implementations is being fueled by the rapid advances and emerging paradigms in HPC at the dawn of exascale computing. Moreover, amidst the explosion of data science, engineering, and medicine, machine learning (ML) integrated with MSM is poised to enhance the capabilities of standard MSM approaches significantly, particularly in the face of increasing problem complexity. The potential to blend MSM, HPC, and ML presents opportunities for unbound innovation and promises to represent the future of MSM and explainable ML that will likely define the fields in the 21st century.

Finding multiple reaction pathways via global optimization of action

Global searching for reaction pathways is a long-standing challenge in computational chemistry and biology. Most existing approaches perform only local searches due to computational complexity. Here we present a computational approach, Action-CSA, to find multiple diverse reaction pathways connecting fixed initial and final states through global optimization of the Onsager-Machlup action using the conformational space annealing (CSA) method. Action-CSA successfully overcomes large energy barriers via crossovers and mutations of pathways and finds all possible pathways of small systems without initial guesses on pathways. The rank order and the transition time distribution of multiple pathways are in good agreement with those of long Langevin dynamics simulations. The lowest action folding pathway of FSD-1 is consistent with recent experiments. The results show that Action-CSA is an efficient and robust computational approach to study the multiple pathways of complex reactions and large-scale conformational changes.

Perspective: Computer simulations of long time dynamics

Atomically detailed computer simulations of complex molecular events attracted the imagination of many researchers in the field as providing comprehensive information on chemical, biological, and physical processes. However, one of the greatest limitations of these simulations is of time scales. The physical time scales accessible to straightforward simulations are too short to address many interesting and important molecular events. In the last decade significant advances were made in different directions (theory, software, and hardware) that significantly expand the capabilities and accuracies of these techniques. This perspective describes and critically examines some of these advances.

Mapping L1 ligase ribozyme conformational switch

L1 ligase (L1L) molecular switch is an in vitro optimized synthetic allosteric ribozyme that catalyzes the regioselective formation of a 5'-to-3' phosphodiester bond, a reaction for which there is no known naturally occurring RNA catalyst. L1L serves as a proof of principle that RNA can catalyze a critical reaction for prebiotic RNA self-replication according to the RNA world hypothesis. L1L crystal structure captures two distinct conformations that differ by a reorientation of one of the stems by around 80Å and are presumed to correspond to the active and inactive state, respectively. It is of great interest to understand the nature of these two states in solution and the pathway for their interconversion. In this study, we use explicit solvent molecular simulation together with a novel enhanced sampling method that utilizes concepts from network theory to map out the conformational transition between active and inactive states of L1L. We find that the overall switching mechanism can be described as a three-state/two-step process. The first step involves a large-amplitude swing that reorients stem C. The second step involves the allosteric activation of the catalytic site through distant contacts with stem C. Using a conformational space network representation of the L1L switch transition, it is shown that the connection between the three states follows different topographical patterns: the stem C swing step passes through a narrow region of the conformational space network, whereas the allosteric activation step covers a much wider region and a more diverse set of pathways through the network.

Experiments and comprehensive simulations of the formation of a helical turn

We investigate the kinetics and thermodynamics of a helical turn formation in the peptide Ac-WAAAH-NH(2). NMR measurements indicate that this peptide has significant tendency to form a structure of a helical turn, while temperature dependent CD establishes the helix fraction at different temperatures. Molecular dynamics and milestoning simulations agree with experimental observables and suggest an atomically detailed picture for the turn formation. Using a network representation, two alternative mechanisms of folding are identified: (i) a direct co-operative mechanism from the unfolded to the folded state without intermediate formation of hydrogen bonds and (ii) an indirect mechanism with structural intermediates with two residues in a helical conformation. This picture is consistent with kinetic measurements that reveal two experimental time scales of sub-nanosecond and several nanoseconds.

AN OVERVIEW OF THE STRETCHED INTERMEDIATE MODEL OF B–Z DNA TRANSITION

Recently, the stretched intermediate model was proposed for the B–Z deoxyribonucleic acid (DNA) transition based on simulation results carried out using the Stochastic Difference Equation (SDE) that showed unwinding and elongation of the oligomer during the transition. This model has proven to be successful in resolving the steric dilemma in short oligomers. However, extending the simulation method to larger DNA strands may prove to be computationally challenging. Such difficulty has led us to adopt a mean field approach using phenomenological interaction potentials to simulate the transition. Like the atomistic approach, the SDE simulations based on the mean field approach, also suggest the presence of a stretched intermediate during the transition.

Revisiting and computing reaction coordinates with Directional Milestoning

The method of Directional Milestoning is revisited. We start from an exact and more general expression and state the conditions and validity of the memory-loss approximation. An algorithm to compute a reaction coordinate from Directional Milestoning data is presented. The reaction coordinate is calculated as a set of discrete jumps between Milestones that maximizes the flux between two stable states. As an application we consider a conformational transition in solvated adenosine. We compare a long molecular dynamic trajectory with Directional Milestoning and discuss the differences between the maximum flux path and minimum energy coordinates.

Onsager-Machlup action-based path sampling and its combination with replica exchange for diffusive and multiple pathways

For sampling multiple pathways in a rugged energy landscape, we propose a novel action-based path sampling method using the Onsager-Machlup action functional. Inspired by the Fourier-path integral simulation of a quantum mechanical system, a path in Cartesian space is transformed into that in Fourier space, and an overdamped Langevin equation is derived for the Fourier components to achieve a canonical ensemble of the path at a finite temperature. To avoid "path trapping" around an initially guessed path, the path sampling method is further combined with a powerful sampling technique, the replica exchange method. The principle and algorithm of our method is numerically demonstrated for a model two-dimensional system with a bifurcated potential landscape. The results are compared with those of conventional transition path sampling and the equilibrium theory, and the error due to path discretization is also discussed.

An experimentally guided umbrella sampling protocol for biomolecules

We present a simple method for utilizing experimental data to improve the efficiency of numerical calculations of free energy profiles from molecular dynamics simulations. The method involves umbrella sampling simulations with restraining potentials based on a known approximate estimate of the free energy profile derived solely from experimental data. The use of the experimental data results in optimal restraining potentials, guides the simulation along relevant pathways, and decreases overall computational time. In demonstration of the method, two systems are showcased. First, guided, unguided (regular) umbrella sampling simulations and exhaustive sampling simulations are compared to each other in the calculation of the free energy profile for the distance between the ends of a pentapeptide. The guided simulation use restraints based on a simulated "experimental" potential of mean force of the end-to-end distance that would be measured by fluorescence resonance energy transfer (obtained from exhaustive sampling). Statistical analysis shows a dramatic improvement in efficiency for a 5 window guided umbrella sampling over 5 and 17 window unguided umbrella sampling simulations. Moreover, the form of the potential of mean force for the guided simulations evolves, as one approaches convergence, along the same milestones as the extensive simulations, but exponentially faster. Second, the method is further validated by replicating the forced unfolding pathway of the titin I27 domain using guiding umbrella sampling potentials determined from actual single molecule pulling data. Comparison with unguided umbrella sampling reveals that the use of guided sampling encourages unfolding simulations to converge faster to a forced unfolding pathway that agrees with previous results and produces a more accurate potential of mean force.

Efficient and verified simulation of a path ensemble for conformational change in a united-residue model of calmodulin

The computational sampling of rare, large-scale, conformational transitions in proteins is a well appreciated challenge-for which a number of potentially efficient path-sampling methodologies have been proposed. Here, we study a large-scale transition in a united-residue model of calmodulin using the "weighted ensemble" (WE) approach of Huber and Kim. Because of the model's relative simplicity, we are able to compare our results with brute-force simulations. The comparison indicates that the WE approach quantitatively reproduces the brute-force results, as assessed by considering (i) the reaction rate, (ii) the distribution of event durations, and (iii) structural distributions describing the heterogeneity of the paths. Importantly, the WE method is readily applied to more chemically accurate models, and by studying a series of lower temperatures, our results suggest that the WE method can increase efficiency by orders of magnitude in more challenging systems.

Exact low-force kinetics from high-force single-molecule unfolding events

Mechanical forces play a key role in crucial cellular processes involving force-bearing biomolecules, as well as in novel single-molecule pulling experiments. We present an exact method that enables one to extrapolate, to low (or zero) forces, entire time-correlation functions and kinetic rate constants from the conformational dynamics either simulated numerically or measured experimentally at a single, relatively higher, external force. The method has twofold relevance: 1), to extrapolate the kinetics at physiological force conditions from molecular dynamics trajectories generated at higher forces that accelerate conformational transitions; and 2), to extrapolate unfolding rates from experimental force-extension single-molecule curves. The theoretical formalism, based on stochastic path integral weights of Langevin trajectories, is presented for the constant-force, constant loading rate, and constant-velocity modes of the pulling experiments. For the first relevance, applications are described for simulating the conformational isomerization of alanine dipeptide; and for the second relevance, the single-molecule pulling of RNA is considered. The ability to assign a weight to each trace in the single-molecule data also suggests a means to quantitatively compare unfolding pathways under different conditions.

A skewed-momenta method to efficiently generate conformational-transition trajectories

We present a novel computational method, the skewed-momenta method (Skew'M), which applies a bias to the Maxwell distribution of initial momenta used to generate ensembles of trajectories. As a result, conformational transitions are accentuated and kinetic properties are calculated more effectively. The connection to the related puddle jumping method is discussed. A reweighting scheme permits the exact calculation of kinetic properties. Applications are presented for the rapid calculation of rate constants for molecular isomerization, and for the efficient reconstruction of free-energy profiles using a straightforward modification of the Jarzynski identity.

Supersymmetric Langevin equation to explore free-energy landscapes

The recently discovered supersymmetric generalizations of the Langevin dynamics and Kramers equation can be utilized for the exploration of free-energy landscapes of systems whose large time-scale separation hampers the usefulness of standard molecular dynamics techniques. The first realistic application is here presented. The system chosen is a minimalist model for a short alanine peptide exhibiting a helix-coil transition.

Transition-event durations in one-dimensional activated processes

Despite their importance in activated processes, transition-event durations--which are much shorter than first passage times--have not received a complete theoretical treatment. The authors therefore study the distribution rhob(t) of durations of transition events over a barrier in a one-dimensional system undergoing overdamped Langevin dynamics. The authors show that rhob(t) is determined by a Fokker-Planck equation with absorbing boundary conditions and obtain a number of results, including (i) the analytic form of the asymptotic short-time transient behavior, which is universal and independent of the potential function; (ii) the first nonuniversal correction to the short-time behavior leading to an estimate of a key physical time scale; (iii) following previous work, a recursive formulation for calculating, exactly, all moments of rhob based solely on the potential function-along with approximations for the distribution based on a small number of moments; and (iv) a high-barrier approximation to the long-time (t-->infinity) behavior of rhob(t). The authors also find that the mean event duration does not depend simply on the barrier-top frequency (curvature) but is sensitive to details of the potential. All of the analytic results are confirmed by transition-path-sampling simulations implemented in a novel way. Finally, the authors discuss which aspects of the duration distribution are expected to be general for more complex systems.

Correct and incorrect nucleotide incorporation pathways in DNA polymerase beta

Tracking the structural and energetic changes in the pathways of DNA replication and repair is central to the understanding of these important processes. Here we report favorable mechanisms of the polymerase-catalyzed phosphoryl transfer reactions corresponding to correct and incorrect nucleotide incorporations in the DNA by using a novel protocol involving energy minimizations, dynamics simulations, quasi-harmonic free energy calculations, and mixed quantum mechanics/molecular mechanics dynamics simulations. Though the pathway proposed may not be unique and invites variations, geometric and energetic arguments support the series of transient intermediates in the phosphoryl transfer pathways uncovered here for both the G:C and G:A systems involving a Grotthuss hopping mechanism of proton transfer between water molecules and the three conserved aspartate residues in pol beta's active-site. In the G:C system, the rate-limiting step is the initial proton hop with a free energy of activation of at least 17 kcal/mol, which corresponds closely to measured k(pol) values. Fidelity discrimination in pol beta can be explained by a significant loss of stability of the closed ternary complex of the enzyme in the G:A system and much higher activation energy of the initial step of nucleophilic attack, namely deprotonation of terminal DNA primer O3'H group. Thus, subtle differences in the enzyme active-site between matched and mismatched base pairs generate significant differences in catalytic performance.

Quantifying the kinetic paths of flexible biomolecular recognition

Biomolecular recognition often involves large conformational changes, sometimes even local unfolding. The identification of kinetic pathways has become a central issue in understanding the nature of binding. A new approach is proposed here to study the dynamics of this binding-folding process through the establishment of a path-integral framework on the underlying energy landscape. The dominant kinetic paths of binding and folding can be determined and quantified. The significant coupling between the binding and folding of biomolecules often exists in many important cellular processes. In this case, the corresponding kinetic paths of binding are shown to be intimately correlated with those of folding and the dynamics becomes quite cooperative. This implies that binding and folding happen concurrently. When the coupling between binding and folding is weak (strong), the kinetic process usually starts with significant folding (binding) first, with the binding (folding) later proceeding to the end. The kinetic rate can be obtained through the contributions from the dominant paths. The rate is shown to have a bell-shaped dependence on temperature in the concentration-saturated regime consistent with experiment. The changes of the kinetics that occur upon changing the parameters of the underlying binding-folding energy landscape are studied.

A line integral reaction path approximation for large systems via nonlinear constrained optimization: Application to alanine dipeptide and the β hairpin of protein G

A variation of the line integral method of Elber with self-avoiding walk has been implemented using a state of the art nonlinear constrained optimization procedure. The new implementation appears to be robust in finding approximate reaction paths for small and large systems. Exact transition states and intermediates for the resulting paths can easily be pinpointed with subsequent application of the conjugate peak refinement method [S. Fischer and M. Karplus, Chem. Phys. Lett. 194, 252 (1992)] and unconstrained minimization, respectively. Unlike previous implementations utilizing a penalty function approach, the present implementation generates an exact solution of the underlying problem. Most importantly, this formulation does not require an initial guess for the path, which makes it particularly useful for studying complex molecular rearrangements. The method has been applied to conformational rearrangements of the alanine dipeptide in the gas phase and in water, and folding of the beta hairpin of protein G in water. In the latter case a procedure was developed to systematically sample the potential energy surface underlying folding and reconstruct folding pathways within the nearest-neighbor hopping approximation.

Dominant kinetic paths on biomolecular binding-folding energy landscape

The identification of kinetic pathways is a central issue in understanding the nature of flexible binding. A new approach is proposed here to study the dynamics of this binding-folding process through the establishment of a path integral framework on the underlying energy landscape. The dominant kinetic paths of binding and folding can be determined and quantified. In this case, the corresponding kinetic paths of binding are shown to be intimately correlated with those of folding and the dynamics becomes quite cooperative. The kinetic time can be obtained through the contributions from the dominant paths and has a U-shape dependence on temperature.

On the calculation of time correlation functions by potential scaling

We present and analyze a general method to calculate time correlation functions from molecular dynamics on scaled potentials for complex systems for which simulation is affected by broken ergodicity. Depending on the value of the scaling factor, correlations can be calculated for times that can be orders of magnitude longer than those accessible to direct simulations. We show that the exact value of the time correlation functions of the original system (i.e., with unscaled potential) can be obtained, in principle, using an action-reweighting scheme based on a stochastic path-integral formalism. Two tests (involving a bistable potential model and a dipeptide bond-vector orientational relaxation) are exemplified to showcase the strengths, as well as the limitations of the approach, and a procedure for the estimation of the time-dependent standard deviation error is outlined.

Energy landscapes and properties of biomolecules

Thermodynamic and dynamic properties of biomolecules can be calculated using a coarse-grained approach based upon sampling stationary points of the underlying potential energy surface. The superposition approximation provides an overall partition function as a sum of contributions from the local minima, and hence functions such as internal energy, entropy, free energy and the heat capacity. To obtain rates we must also sample transition states that link the local minima, and the discrete path sampling method provides a systematic means to achieve this goal. A coarse-grained picture is also helpful in locating the global minimum using the basin-hopping approach. Here we can exploit a fictitious dynamics between the basins of attraction of local minima, since the objective is to find the lowest minimum, rather than to reproduce the thermodynamics or dynamics.

Quantifying kinetic paths of protein folding

We propose a new approach to activated protein folding dynamics via a diffusive path integral framework. The important issues of kinetic paths in this situation can be directly addressed. This leads to the identification of the kinetic paths of the activated folding process, and provides a direct tool and language for the theoretical and experimental community to understand the problem better. The kinetic paths giving the dominant contributions to the long-time folding activation dynamics can be quantitatively determined. These are shown to be the instanton paths. The contributions of these instanton paths to the kinetics lead to the "bell-like" shape folding rate dependence on temperature, which is in good agreement with folding kinetic experiments and simulations. The connections to other approaches as well as the experiments of the protein folding kinetics are discussed.

The stretched intermediate model of B-Z DNA transition

There have been numerous attempts to describe the mechanism of B-Z transition. Our simulations based on the stochastic difference equation with length algorithm show that a short DNA oligomer will tend to unwind and overstretch during the transition. The overstretching of DNA is then understood from the Zhou, Zhang, and Ou-Yang model. Unlike the Harvey model, the stretched intermediate model does not pose any steric dilemma; we are able to show that the chain sense reversal progresses spontaneously using the stretched intermediate model. A nonlinear DNA model is used to describe the origins and mechanism of base rotation in the stretched intermediate state of DNA. We also propose an experiment that can verify the existence of a stretched intermediate state during B-Z transition, thus opening up fresh grounds for experimentation in this field.

Atomically detailed simulations of helix formation with the stochastic difference equation

An algorithm is described to compute approximate classical trajectories as a boundary value problem with an integration step in the arc length. High-frequency motions are filtered out when a large integration step is used, maintaining the stability of the algorithm. At the limit of high filtering (large steps), but still offering an accurate description of the continuous path, the trajectory approaches the steepest descent path (SDP). The steepest descent path is widely used as a reaction coordinate in chemical systems. At intermediate step sizes, some inertial motions remain, interpolating between reaction coordinates and exact classical trajectories. Numerical studies of spatial and energetic properties of meta-trajectories are carried out. Two systems are considered here: valine dipeptide and the folding of a small helical protein. Although thermodynamic properties of meta-trajectories are affected by the filtering, the ordering of events remains similar for substantial differences in trajectory resolution.

Kinetics of cytochrome C folding: atomically detailed simulations

The vast range of time scales (from nanoseconds to seconds) during protein folding is a challenge for experiments and computations. To make concrete predictions on folding mechanisms, atomically detailed simulations of protein folding, using potentials derived from chemical physics principles, are desired. However, due to their computational complexity, straightforward molecular dynamics simulations of protein folding are impossible today. An alternative algorithm is used that makes it possible to compute approximate atomically detailed long time trajectories (the Stochastic Difference Equation in Length). This algorithm is used to compute 26 atomically detailed folding trajectories of cytochrome c (a millisecond process). The early collapse of the protein chain (with marginal formation of secondary structure), and the earlier formation of the N and C helices (compare to the 60's helix) are consistent with the experiment. The existence of an energy barrier upon entry to the molten globule is examined as well. In addition to (favorable) comparison to experiments, we show that non-native contacts drive the formation of the molten globule. In contrast to popular folding models, the non-native contacts do not form off-pathway kinetic traps in cytochrome c.

Thermal unfolding simulations of apo-calmodulin using leap-dynamics

The simulation method leap-dynamics (LD) has been applied to protein thermal unfolding simulations to investigate domain-specific unfolding behavior. Thermal unfolding simulations of the 148-residue protein apo-calmodulin with implicit solvent were performed at temperatures 290 K, 325 K, and 360 K and compared with the corresponding molecular dynamics trajectories in terms of a number of calculated conformational parameters. The main experimental results of unfolding are reproduced in showing the lower stability of the C-domain: at 290 K, both the N- and C-domains are essentially stable; at 325 K, the C-domain unfolds, whereas the N-domain remains folded; and at 360 K, both domains unfold extensively. This behavior could not be reproduced by molecular dynamics simulations alone under the same conditions. These results show an encouraging degree of convergence between experiment and LD simulation. The simulations are able to describe the overall plasticity of the apo-calmodulin structure and to reveal details such as reversible folding/unfolding events within single helices. The results show that by using the combined application of a fast and efficient sampling routine with a detailed molecular dynamics force field, unfolding simulations of proteins at atomic resolution are within the scope of current computational power.

Ion permeation through the gramicidin channel: atomically detailed modeling by the Stochastic Difference Equation

Atomically detailed descriptions of ionic solution, membrane, and the gramicidin channel are used to compute molecular dynamics trajectories of ion permeation. The microsecond trajectories are calculated with the Stochastic Difference Equation (SDE), which provides approximate solutions to the equations of motions (with filtered high-frequency modes) of extended timescales. The relative permeations of lithium, sodium, and potassium are estimated by using a novel, kinetic cycle protocol and are compared with experiment. The transport through native gramicidin and one fluoro-valine variant is considered as well. Qualitative agreement between theory and experiment is obtained. The faster permeation rate of sodium compared to lithium is reproduced in the calculations. The calculations also reproduce the slower diffusion through a gramicidin with fluorinated valine compared to native gramicidin. The calculations are inconclusive about the relative rates of potassium and sodium. The experiment suggests that potassium permeates more quickly. We directly probe the kinetics of a biophysical process at a relevant time window without reducing the atomically detailed description of the system. The calculations were able to capture subtle balances between binding and diffusion that determine permeation rates. The same model gave the correct ordering of diffusion rates for cases in which electrostatic binding has opposite effects and must be supplemented by dynamic factors. Diffusion rates are faster when favorable electrostatic interactions of ions in the channel (compared to the solvent) are observed. Studies of a gramicidin variant suggest an opposite effect, in which permeation is faster for the less polar channel, indicating dynamic effects. Although both trends can be explained qualitatively, it is not possible to predict (before doing the SDE calculations) which factor is more important.

Parallel-in-time molecular-dynamics simulations

While there have been many progress in the field of multiscale simulations in the space domain, in particular, due to efficient parallelization techniques, much less is known in the way to perform similar approaches in the time domain. In this paper we show on two examples that, provided we can describe in a rough but still accurate way the system under consideration, it is indeed possible to parallelize molecular dynamics simulations in time by using the recently introduced pararealalgorithm. The technique is most useful for ab initio simulations.

Long time dynamics of complex systems

Molecular dynamics trajectories of large biological molecules are restricted to nanoseconds. We describe a computational method, based on optimization of a functional, to extend the time of molecular simulations by orders of magnitude. Variants of our technique have already produced microsecond and millisecond trajectories. The large steps enable feasible computations of atomically detailed approximate trajectories. Numerical examples are provided: (i) a conformational change in blocked glycine peptide and (ii) helix formation of an alanine-rich peptide.