Quantifying kinetic paths of protein folding

Perspectives on the landscape and flux theory for describing emergent behaviors of the biological systems

We give a review on the landscape theory of the equilibrium biological systems and landscape-flux theory of the nonequilibrium biological systems as the global driving force. The emergences of the behaviors, the associated thermodynamics in terms of the entropy and free energy and dynamics in terms of the rate and paths have been quantitatively demonstrated. The hierarchical organization structures have been discussed. The biological applications ranging from protein folding, biomolecular recognition, specificity, biomolecular evolution and design for equilibrium systems as well as cell cycle, differentiation and development, cancer, neural networks and brain function, and evolution for nonequilibrium systems, cross-scale studies of genome structural dynamics and experimental quantifications/verifications of the landscape and flux are illustrated. Together, this gives an overall global physical and quantitative picture in terms of the landscape and flux for the behaviors, dynamics and functions of biological systems.

Noise-induced quasiperiod and period switching

We employ a typical genetic circuit model to explore how noise can influence dynamic structure. With the increase of a key interactive parameter, the model will deterministically go through two bifurcations and three dynamic structure regions. We find that a quasiperiodic component, which is not allowed by deterministic dynamics, will be generated by noise inducing in the first two regions, and this quasiperiod will be more and more stable along with the increase in noise. In particular, in the second region the quasiperiod will compete with a stable limit cycle and perform a new transient rhythm. Furthermore, we ascertain the entropy production rate and the heat dissipation rate, and discover a minimal value with theoretical elucidation. In the end, we unveil the mechanism of the formation of quasiperiods, and show a practical biological example. We expect this work to be helpful in solving some biological or ecological problems, such as the genetic origin of periodical cicadas and population dynamics with fluctuation.

Physics of biomolecular recognition and conformational dynamics

Biomolecular recognition usually leads to the formation of binding complexes, often accompanied by large-scale conformational changes. This process is fundamental to biological functions at the molecular and cellular levels. Uncovering the physical mechanisms of biomolecular recognition and quantifying the key biomolecular interactions are vital to understand these functions. The recently developed energy landscape theory has been successful in quantifying recognition processes and revealing the underlying mechanisms. Recent studies have shown that in addition to affinity, specificity is also crucial for biomolecular recognition. The proposed physical concept of intrinsic specificity based on the underlying energy landscape theory provides a practical way to quantify the specificity. Optimization of affinity and specificity can be adopted as a principle to guide the evolution and design of molecular recognition. This approach can also be used in practice for drug discovery using multidimensional screening to identify lead compounds. The energy landscape topography of molecular recognition is important for revealing the underlying flexible binding or binding-folding mechanisms. In this review, we first introduce the energy landscape theory for molecular recognition and then address four critical issues related to biomolecular recognition and conformational dynamics: (1) specificity quantification of molecular recognition; (2) evolution and design in molecular recognition; (3) flexible molecular recognition; (4) chromosome structural dynamics. The results described here and the discussions of the insights gained from the energy landscape topography can provide valuable guidance for further computational and experimental investigations of biomolecular recognition and conformational dynamics.

Chemical Langevin equation: A path-integral view of Gillespie's derivation

In 2000, Gillespie rehabilitated the chemical Langevin equation (CLE) by describing two conditions that must be satisfied for it to yield a valid approximation of the chemical master equation (CME). In this work, we construct an original path-integral description of the CME and show how applying Gillespie's two conditions to it directly leads to a path-integral equivalent to the CLE. We compare this approach to the path-integral equivalent of a large system size derivation and show that they are qualitatively different. In particular, both approaches involve converting many sums into many integrals, and the difference between the two methods is essentially the difference between using the Euler-Maclaurin formula and using Riemann sums. Our results shed light on how path integrals can be used to conceptualize coarse-graining biochemical systems and are readily generalizable.

Quantifying the landscape and kinetic paths for epithelial-mesenchymal transition from a core circuit

Epithelial-mesenchymal transition (EMT), as a crucial process in embryonic development and cancer metastasis, has been investigated extensively. However, how to quantify the global stability and transition dynamics for EMT under fluctuations remains to be elucidated. Starting from a core EMT genetic circuit composed of three key proteins or microRNAs (microRNA-200, ZEB and SNAIL), we uncovered the potential landscape for the EMT process. Three attractors emerge from the landscape, which correspond to epithelial, mesenchymal and partial EMT states respectively. Based on the landscape, we analyzed two important quantities of the EMT system: the barrier heights between different basins of attraction that describe the degree of difficulty for EMT or backward transition, and the mean first passage time (MFPT) that characterizes the kinetic transition rate. These quantities can be harnessed as measurements for the stability of cell types and the degree of difficulty of transitions between different cell types. We also calculated the minimum action paths (MAPs) by path integral approaches. The MAP delineates the transition processes between different cell types quantitatively. We propose two different EMT processes: a direct EMT from E to P, and a step-wise EMT going through an intermediate state, depending on different extracellular environments. The landscape and kinetic paths we acquired offer a new physical and quantitative way for understanding the mechanisms of EMT processes, and indicate the possible roles for the intermediate states.

Exploring multi-dimensional coordinate-dependent diffusion dynamics on the energy landscape of protein conformation change

We explore the multi-dimensional diffusion dynamics of protein conformational change. We found in general that the diffusion is anisotropic and inhomogeneous. The directional and positional dependence of diffusion have significant impacts on the protein conformational kinetics: the dominant kinetic path of conformational change is shifted from the naively expected steepest decent gradient paths. The kinetic transition state is shifted away from the transition state. The effective kinetic free energy barrier height, determining the kinetic rate of the conformational change, is shifted away from the one estimated from the thermodynamic free energy barrier. The shift of the transition state in position and value will modify the phi value analysis for identification of hot residues and interactions responsible for conformational dynamics. Ongoing and future experiments can test the predictions of the model.

Non-equilibrium transition state rate theory

Transition state or Kramers' rate theory has been used to quantify the kinetic speed of many chemical, physical and biological equilibrium processes successfully.

The dominant folding route minimizes backbone distortion in SH3

Energetic frustration in protein folding is minimized by evolution to create a smooth and robust energy landscape. As a result the geometry of the native structure provides key constraints that shape protein folding mechanisms. Chain connectivity in particular has been identified as an essential component for realistic behavior of protein folding models. We study the quantitative balance of energetic and geometrical influences on the folding of SH3 in a structure-based model with minimal energetic frustration. A decomposition of the two-dimensional free energy landscape for the folding reaction into relevant energy and entropy contributions reveals that the entropy of the chain is not responsible for the folding mechanism. Instead the preferred folding route through the transition state arises from a cooperative energetic effect. Off-pathway structures are penalized by excess distortion in local backbone configurations and contact pair distances. This energy cost is a new ingredient in the malleable balance of interactions that controls the choice of routes during protein folding.

Fluctuations in the ensemble of reaction pathways

The dominant reaction pathway is a rigorous framework to microscopically compute the most probable trajectories, in nonequilibrium transitions. In the low-temperature regime, such dominant pathways encode the information about the reaction mechanism and can be used to estimate nonequilibrium averages of arbitrary observables. On the other hand, at sufficiently high temperatures, the stochastic fluctuations around the dominant paths become important and have to be taken into account. In this work, we develop a technique to systematically include the effects of such stochastic fluctuations, to order k(B)T. This method is used to compute the probability for a transition to take place through a specific reaction channel and to evaluate the reaction rate.

Multi-scaled explorations of binding-induced folding of intrinsically disordered protein inhibitor IA3 to its target enzyme

Biomolecular function is realized by recognition, and increasing evidence shows that recognition is determined not only by structure but also by flexibility and dynamics. We explored a biomolecular recognition process that involves a major conformational change – protein folding. In particular, we explore the binding-induced folding of IA3, an intrinsically disordered protein that blocks the active site cleft of the yeast aspartic proteinase saccharopepsin (YPrA) by folding its own N-terminal residues into an amphipathic alpha helix. We developed a multi-scaled approach that explores the underlying mechanism by combining structure-based molecular dynamics simulations at the residue level with a stochastic path method at the atomic level. Both the free energy profile and the associated kinetic paths reveal a common scheme whereby IA3 binds to its target enzyme prior to folding itself into a helix. This theoretical result is consistent with recent time-resolved experiments. Furthermore, exploration of the detailed trajectories reveals the important roles of non-native interactions in the initial binding that occurs prior to IA3 folding. In contrast to the common view that non-native interactions contribute only to the roughness of landscapes and impede binding, the non-native interactions here facilitate binding by reducing significantly the entropic search space in the landscape. The information gained from multi-scaled simulations of the folding of this intrinsically disordered protein in the presence of its binding target may prove useful in the design of novel inhibitors of aspartic proteinases.

The intrinsically disordered peptide IA3 is the endogenous inhibitor for the enzyme named yeast aspartic proteinase saccharopepsin (YPrA). In the presence of YPrA, IA3 folds itself into an amphipathic helix that blocks the active site cleft of the enzyme. We developed a multi-scaled approach to explore the underlying mechanism of this binding-induced ordering transition. Our approach combines a structure-based molecular dynamics model at the residue level with a stochastic path method at the atomic level. Our simulations suggest that IA3 inhibits YPrA through an induced-fit mechanism where the enzyme (YPrA) induces conformational change of its inhibitor (IA3). This expands the definition of an induced-fit model from its original meaning that the binding of substrate (IA3) drives conformational change in the protein (YPrA). Our result is consistent with recent kinetic experiments and provides a microscopic explanation for the underlying mechanism. We also discuss the important roles of non-native interactions and backtracking. These results enrich our understanding of the enzyme-inhibition mechanism and may have value in the design of drugs.

Protein functional landscapes, dynamics, allostery: a tortuous path towards a universal theoretical framework

Energy landscape theories have provided a common ground for understanding the protein folding problem, which once seemed to be overwhelmingly complicated. At the same time, the native state was found to be an ensemble of interconverting states with frustration playing a more important role compared to the folding problem. The landscape of the folded protein – the native landscape – is glassier than the folding landscape; hence, a general description analogous to the folding theories is difficult to achieve. On the other hand, the native basin phase volume is much smaller, allowing a protein to fully sample its native energy landscape on the biological timescales. Current computational resources may also be used to perform this sampling for smaller proteins, to build a ‘topographical map’ of the native landscape that can be used for subsequent analysis. Several major approaches to representing this topographical map are highlighted in this review, including the construction of kinetic networks, hierarchical trees and free energy surfaces with subsequent structural and kinetic analyses. In this review, we extensively discuss the important question of choosing proper collective coordinates characterizing functional motions. In many cases, the substates on the native energy landscape, which represent different functional states, can be used to obtain variables that are well suited for building free energy surfaces and analyzing the protein's functional dynamics. Normal mode analysis can provide such variables in cases where functional motions are dictated by the molecule's architecture. Principal component analysis is a more expensive way of inferring the essential variables from the protein's motions, one that requires a long molecular dynamics simulation. Finally, the two popular models for the allosteric switching mechanism, ‘preexisting equilibrium’ and ‘induced fit’, are interpreted within the energy landscape paradigm as extreme points of a continuum of transition mechanisms. Some experimental evidence illustrating each of these two models, as well as intermediate mechanisms, is presented and discussed.

Minimal folding pathways for coarse-grained biopolymer fragments

The minimal folding pathway or trajectory for a biopolymer can be defined as the transformation that minimizes the total distance traveled between a folded and an unfolded structure. This involves generalizing the usual Euclidean distance from points to one-dimensional objects such as a polymer. We apply this distance here to find minimal folding pathways for several candidate protein fragments, including the helix, the beta-hairpin, and a nonplanar structure where chain noncrossing is important. Comparing the distances traveled with root mean-squared distance and mean root-squared distance, we show that chain noncrossing can have large effects on the kinetic proximity of apparently similar conformations. Structures that are aligned to the beta-hairpin by minimizing mean root-squared distance, a quantity that closely approximates the true distance for long chains, show globally different orientation than structures aligned by minimizing root mean-squared distance.

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.

Folding time predictions from all-atom replica exchange simulations

We present an approach to predicting the folding time distribution from all-atom replica exchange simulations. This is accomplished by approximating the multidimensional folding process as stochastic reaction-coordinate dynamics for which effective drift velocities and diffusion coefficients are determined from the short-time replica exchange simulations. Our approach is applied to the folding of the second beta-hairpin of the B domain of protein G. The folding time prediction agrees quite well with experimental measurements. Therefore, we have in hand a fast numerical tool for calculating the folding kinetic properties from all-atom "first-principles" models.

MinActionPath: maximum likelihood trajectory for large-scale structural transitions in a coarse-grained locally harmonic energy landscape

The non-linear problem of simulating the structural transition between two known forms of a macromolecule still remains a challenge in structural biology. The problem is usually addressed in an approximate way using 'morphing' techniques, which are linear interpolations of either the Cartesian or the internal coordinates between the initial and end states, followed by energy minimization. Here we describe a web tool that implements a new method to calculate the most probable trajectory that is exact for harmonic potentials; as an illustration of the method, the classical Calpha-based Elastic Network Model (ENM) is used both for the initial and the final states but other variants of the ENM are also possible. The Langevin equation under this potential is solved analytically using the Onsager and Machlup action minimization formalism on each side of the transition, thus replacing the original non-linear problem by a pair of linear differential equations joined by a non-linear boundary matching condition. The crossover between the two multidimensional energy curves around each state is found numerically using an iterative approach, producing the most probable trajectory and fully characterizing the transition state and its energy. Jobs calculating such trajectories can be submitted on-line at: http://lorentz.dynstr.pasteur.fr/joel/index.php.

Temperature-dependent folding pathways of Pin1 WW domain: an all-atom molecular dynamics simulation of a Gō model

We study the folding thermodynamics and kinetics of the Pin1 WW domain, a three-stranded beta-sheet protein, by using all-atom (except nonpolar hydrogens) discontinuous molecular dynamics simulations at various temperatures with a Gō model. The protein exhibits a two-state folding kinetics near the folding transition temperature. A good agreement between our simulations and the experimental measurements by the Gruebele group has been found, and the simulation sheds new insights into the structure of transition state, which is hard to be straightforwardly captured in experiments. The simulation also reveals that the folding pathways at approximately the transition temperature and at low temperatures are much different, and an intermediate state at a low temperature is predicted. The transition state of this small beta-protein at its folding transition temperature has a well-established hairpin 1 made of beta1 and beta2 strands while its low-temperature kinetic intermediate has a formed hairpin 2 composed of beta2 and beta3 strands. Theoretical results are compared with other simulation results as well as available experimental data. This study confirms that specific side-chain packing in an all-atom Gō model can yield a reasonable prediction of specific folding kinetics for a given protein. Different folding behaviors at different temperatures are interpreted in terms of the interplay of entropy and enthalpy in folding process.

Dominant pathways in protein folding

We present a method to investigate the kinetics of protein folding and the dynamics underlying the formation of secondary and tertiary structures during the entire reaction. By writing the solution of the Fokker-Planck equation in terms of a path integral, we derive a Hamilton-Jacobi variational principle from which we are able to compute the most probable pathway of folding. The method is applied to the folding of the Villin headpiece subdomain simulated using a Go model. An initial collapsing phase driven by the initial configuration is followed by a rearrangement phase, in which secondary structures are formed and all computed paths display strong similarities. This completely general method does not require the prior knowledge of any reaction coordinate and is an efficient tool to perform simulations of the entire folding process with available computers.

Effective stochastic dynamics on a protein folding energy landscape

We present an approach to protein folding kinetics using stochastic reaction-coordinate dynamics, in which the effective drift velocities and diffusion coefficients are determined from microscopic simulation data. The resultant Langevin equation can then be used to directly simulate the folding process. Here, we test this approach by applying it to a toy two-state dynamical system and to a funnellike structure-based (Go-type) model. The folding time predictions agree very well with full simulation results. Therefore, we have in hand a fast numerical tool for calculating the folding kinetic properties, even when full simulations are not feasible. In addition, the local drift and diffusion coefficients provide an alternative way to compute the free energy profile in cases where only local sampling can be achieved.

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.

Kinetic definition of protein folding transition state ensembles and reaction coordinates

Using distributed molecular dynamics simulations we located four distinct folding transitions for a 39-residue betabetaalphabeta protein fold. To characterize the nature of each room temperature transition, we calculated the probability of transmission for 500 points along each free energy barrier. We introduced a method for determining transition states by employing the transmission probability, Ptrans, and determined which conformations were transition state ensemble members (Ptrans approximately 0.5). The transmission probability may be used to characterize the barrier in several ways. For example, we ran simulations at 82 degrees C, determined the change in Ptrans with temperature for all 2,000 conformations, and quantified Hammond behavior directly using Ptrans correlation. Additionally, we propose that diffusion along Ptrans may provide the configurational diffusion rate at the top of the barrier. Specifically, given a transition state conformation x0 with estimated Ptrans=0.5, we selected a large set of subsequent conformations from independent trajectories, each exactly a small time deltat after x0 (250 ps). Calculating Ptrans for the new trial conformations, we generated the P(Ptrans|deltat=250 ps) distribution that reflected diffusion. This approach provides a novel perspective on the diffusive nature of a protein folding transition and provides a framework for a quantitative study of activated relaxation kinetics.