Reaction coordinates and rates from transition paths

Is protein folding sub-diffusive?

Protein folding dynamics is often described as diffusion on a free energy surface considered as a function of one or few reaction coordinates. However, a growing number of experiments and models show that, when projected onto a reaction coordinate, protein dynamics is sub-diffusive. This raises the question as to whether the conventionally used diffusive description of the dynamics is adequate. Here, we numerically construct the optimum reaction coordinate for a long equilibrium folding trajectory of a Go model of a -repressor protein. The trajectory projected onto this coordinate exhibits diffusive dynamics, while the dynamics of the same trajectory projected onto a sub-optimal reaction coordinate is sub-diffusive. We show that the higher the (cut-based) free energy profile for the putative reaction coordinate, the more diffusive the dynamics become when projected on this coordinate. The results suggest that whether the projected dynamics is diffusive or sub-diffusive depends on the chosen reaction coordinate. Protein folding can be described as diffusion on the free energy surface as function of the optimum reaction coordinate. And conversely, the conventional reaction coordinates, even though they might be based on physical intuition, are often sub-optimal and, hence, show sub-diffusive dynamics.

To understand dynamics of complex systems with many degrees of freedom, one often projects it onto one or several collective variables. Protein folding, the complex, concerted motion of a protein chain towards a unique three-dimensional structure, is one example of where such reduction of complexity is useful. It is usually assumed that the projected dynamics is diffusive. However, many experiments and simulations have shown that the projected dynamics is sub-diffusive, i.e., the mean square displacement grows slower than linear with time. It means that the dynamics has a memory; that the free energy surface together with diffusion coefficient do not properly define the dynamics; and that such projections cannot be used to accurately describe dynamics. Here, we show that if one carefully constructs the reaction coordinate by optimizing (maximizing) its free energy profile, one can use a simple (memory-less) diffusive description. Loosely speaking, when the complex dynamics is projected onto a simple coordinate, all the complexity of the original dynamics goes into the memory of the projected dynamics. If the dynamics is projected onto the (complex) optimum reaction coordinate, all the complexity of the original dynamics is in the reaction coordinate, and the projected dynamics is simple.

Extracting physically intuitive reaction coordinates from transition networks of a beta-sheet miniprotein

Simulations are important for understanding complex reactions, but their interpretation is challenging owing to the large number of degrees of freedom typically involved. To address this issue, various means for relating the dynamics of a stochastic system to its structural and energetic features have been introduced. Here, we show how two leading approaches can be combined to advantage. We use the network of transitions observed in a reversible folding/unfolding simulation of a 20-residue three-stranded antiparallel beta-sheet peptide (beta3s) to estimate the probabilities of committing to stable states (the native state and major nonnative states), and these then serve as the basis for an efficient statistical procedure for identifying physical variables that describe the dynamics. We find that a single coordinate that jointly characterizes the formation of the two native turns of beta3s can adequately describe the overall folding process, despite its complex nature. Additional features associated with major pathways leading from individual nonnative states are resolved; indeed, a key result is an improved understanding of the unfolded state. Connections to other methods for analyzing complex reactions are discussed.

Nonequilibrium solvent effects in reaction kinetics — Steady-state solutions for the Agmon–Hopfield two-dimensional stochastic model

We propose an approximate localized-sink approach to the solution of the Agmon–Hopfield two-dimensional stochastic model for reactions in viscous media. The approach yields simple expressions for the sink location and the reaction rate constant, easy to use in the analysis of experimental data, and allows an intuitive phenomenological interpretation of the overall process in terms of a two-step kinetic scheme.

Systematic determination of order parameters for chain dynamics using diffusion maps

We employ the diffusion map approach as a nonlinear dimensionality reduction technique to extract a dynamically relevant, low-dimensional description of n-alkane chains in the ideal-gas phase and in aqueous solution. In the case of C8 we find the dynamics to be governed by torsional motions. For C16 and C24 we extract three global order parameters with which we characterize the fundamental dynamics, and determine that the low free-energy pathway of globular collapse proceeds by a "kink and slide" mechanism, whereby a bend near the end of the linear chain migrates toward the middle to form a hairpin and, ultimately, a coiled helix. The low-dimensional representation is subtly perturbed in the solvated phase relative to the ideal gas, and its geometric structure is conserved between C16 and C24. The methodology is directly extensible to biomolecular self-assembly processes, such as protein folding.

Cooperativity at different space and time scales in multiscale protein dynamics

A method proposed by Matsunaga [Phys. Rev. Lett. 99, 238103 (2007)] is applied to simple stochastic models and two model proteins composed of 46 amino beads with three different kinds of residues. The method, which is based on the combination of the principal component analysis and the finite size Lyapunov exponent, characterize the coarse-grained dynamics in different spatiotemporal hierarchies in protein dynamics. The application of the method to model proteins reveals that the low-indexed (large-variance) principal components carry less-divergent, regularized dynamics at the coarse-grained scales on a less-frustrated energy landscape, whereas this less-divergent nature is less pronounced for a protein model with a more frustrated energy landscape. It is also revealed that our technique can differentiate the collective motions on the projected principal component space inherent to the system and the apparent collective behavior which can appear even in high-dimensional stochastic systems.

Dependence of protein folding stability and dynamics on the density and composition of macromolecular crowders

We investigate the effect of macromolecular crowding on protein folding, using purely repulsive crowding particles and a self-organizing polymer model of protein folding. We find that the variation in folding stability with crowder size for typical alpha-, beta-, and alpha/beta-proteins is well described by an adaptation of the scaled particle theory. The native state, the transition state, and the unfolded protein are treated as effective hard spheres, with the folded and transition state radii independent of the size and concentration of the crowders. Remarkably, we find that, as the effective unfolded state radius is very weakly dependent on the crowder concentration, it can also be approximated by a single size. The same model predicts the effect of crowding on the folding barrier and therefore refolding rates with no adjustable parameters. A simple extension of the scaled-particle theory model, assuming additivity, can also describe the behavior of mixtures of crowding particles.

Many local motions cooperate to produce the adenylate kinase conformational transition

Conformational transitions are functionally important in many proteins. In the enzyme adenylate kinase (AK), two small domains (LID and NMP) close over the larger CORE domain; the reverse (opening) motion limits the rate of catalytic turnover. Here, using double-well Gō simulations of Escherichia coli AK, we elaborate on previous investigations of the AK transition mechanism by characterizing the contributions of rigid-body (Cartesian), backbone dihedral, and contact motions to transition-state (TS) properties. In addition, we compare an apo simulation to a pseudo-ligand-bound simulation to reveal insights into allostery. In Cartesian space, LID closure precedes NMP closure in the bound simulation, consistent with prior coarse-grained models of the AK transition. However, NMP-first closure is preferred in the apo simulation. In backbone dihedral space, we find that, as expected, backbone fluctuations are reduced in the O/C transition in parts of all three domains. Among these "quenching" residues, most in the CORE, especially residues 11-13, are rigidified in the TS of the bound simulation, while residues 42-44 in the NMP are flexible in the TS. In contact space, in both apo and bound simulations, one nucleus of closed-state contacts includes parts of the NMP and CORE; CORE-LID contacts are absent in the TS of the apo simulation but formed in the TS of the bound simulation. From these results, we predict mutations that will perturb the opening and/or closing transition rates by changing the entropy of dihedrals and/or the enthalpy of contacts. Furthermore, regarding allostery, the fully closed structure is populated in the apo simulation, but our contact results imply that ligand binding shifts the preferred O/C transition pathway, thus precluding a simple conformational selection mechanism. Finally, the analytical approach and the insights derived from this work may inform the rational design of flexibility and allostery in proteins.

Predicting the reaction coordinates of millisecond light-induced conformational changes in photoactive yellow protein

Understanding the dynamics of large-scale conformational changes in proteins still poses a challenge for molecular simulations. We employ transition path sampling of explicit solvent molecular dynamics trajectories to obtain atomistic insight in the reaction network of the millisecond timescale partial unfolding transition in the photocycle of the bacterial sensor photoactive yellow protein. Likelihood maximization analysis predicts the best model for the reaction coordinates of each substep as well as tentative transition states, without further simulation. We find that the unfolding of the alpha-helical region 43-51 is followed by sequential solvent exposure of both Glu46 and the chromophore. Which of these two residues is exposed first is correlated with the presence of a salt bridge that is part of the N-terminal domain. Additional molecular dynamics simulations indicate that the exposure of the chromophore does not result in a productive pathway. We discuss several possibilities for experimental validation of these predictions. Our results open the way for studying millisecond conformational changes in other medium-sized (signaling) proteins.

Exploration of effective potential landscapes using coarse reverse integration

We describe a reverse integration approach for the exploration of low-dimensional effective potential landscapes. Coarse reverse integration initialized on a ring of coarse states enables efficient navigation on the landscape terrain: Escape from local effective potential wells, detection of saddle points, and identification of significant transition paths between wells. We consider several distinct ring evolution modes: Backward stepping in time, solution arc length, and effective potential. The performance of these approaches is illustrated for a deterministic problem where the energy landscape is known explicitly. Reverse ring integration is then applied to noisy problems where the ring integration routine serves as an outer wrapper around a forward-in-time inner simulator. Two versions of such inner simulators are considered: A Gillespie-type stochastic simulator and a molecular dynamics simulator. In these "equation-free" computational illustrations, estimation techniques are applied to the results of short bursts of inner simulation to obtain the unavailable (in closed-form) quantities (local drift and diffusion coefficient estimates) required for reverse ring integration; this naturally leads to approximations of the effective landscape.

Coordinate-dependent diffusion in protein folding

Diffusion on a low-dimensional free-energy surface is a remarkably successful model for the folding dynamics of small single-domain proteins. Complicating the interpretation of both simulations and experiments is the expectation that the effective diffusion coefficient D will in general depend on the position along the folding coordinate, and this dependence may vary for different coordinates. Here we explore the position dependence of D, its connection to protein internal friction, and the consequences for the interpretation of single-molecule experiments. We find a large decrease in D from unfolded to folded, for reaction coordinates that directly measure fluctuations in Cartesian configuration space, including those probed in single-molecule experiments. In contrast, D is almost independent of Q, the fraction of native amino acid contacts: Near the folded state, small fluctuations in position cause large fluctuations in Q, and vice versa for the unfolded state. In general, position-dependent free energies and diffusion coefficients for any two good reaction coordinates that separate reactant, product, and transition states, are related by a simple transformation, as we demonstrate. With this transformation, we obtain reaction coordinates with position-invariant D. The corresponding free-energy surfaces allow us to justify the assumptions used in estimating the speed limit for protein folding from experimental measurements of the reconfiguration time in the unfolded state, and also reveal intermediates hidden in the original free-energy projection. Lastly, we comment on the design of future single-molecule experiments that probe the position dependence of D directly.

Forward flux sampling for rare event simulations

Rare events are ubiquitous in many different fields, yet they are notoriously difficult to simulate because few, if any, events are observed in a conventional simulation run. Over the past several decades, specialized simulation methods have been developed to overcome this problem. We review one recently developed class of such methods, known as forward flux sampling. Forward flux sampling uses a series of interfaces between the initial and final states to calculate rate constants and generate transition paths for rare events in equilibrium or nonequilibrium systems with stochastic dynamics. This review draws together a number of recent advances, summarizes several applications of the method and highlights challenges that remain to be overcome.

Nearly symmetrical proteins: folding pathways and transition states

The folding pathways of the B domain of protein A have been the subject of many experimental and computational studies. Based on a statistical mechanical model, it has been suggested that the native state symmetry leads to multiple pathways, highly dependent on temperature and denaturant concentration. Experiments, however, have not confirmed this scenario. By considering four nearly symmetrical proteins, one of them being the above molecule, here we show that, if contact energies are properly taken into account, a different picture emerges from kinetic simulations of the above-mentioned model. This is characterized by a dominant folding pathway, which is consistent with the most recent experimental results. Given the simplicity of the model, we also report on a direct sampling of the transition state.

Orientational dynamics and dielectric response of nanopore water

We present numerical calculations, simulation results, and analytical considerations for the frequency-dependent dielectric constant of single-file water in narrow nanopores, described by a recently developed dipole lattice model. We find Debye relaxation over all length scales with relaxation times that strongly depend on pore length. This behavior is analyzed in terms of the dynamics of orientational defects leading to simple quantitative expressions for the static dielectric susceptibility and the relaxation time in the limits of short and long pores. Based on these formulas, we suggest how the predicted macroscopic order of nanopore water can be probed via dielectric spectroscopy and explain how the excitation energy, diffusion constant, and effective interaction of the defects that destroy the order can be extracted from such measurements.

Transition path sampling and forward flux sampling. Applications to biological systems

The last decade has seen a rapid growth in the number of simulation methods and applications dealing with the sampling of transition pathways of rare nanoscale events. Such studies are crucial, for example, for understanding the mechanism and kinetics of conformational transitions and enzymatic events associated with the function of biomolecules. In this review, a broad account of transition path sampling approaches is provided, starting from the general concepts, progressing to the specific principles that underlie some of the most important methods, and eventually singling out the so-called forward flux sampling method for a more detailed description. This is done because forward flux sampling, despite its appealing simplicity and potential efficiency, has thus far received limited attention from practitioners. While path sampling methods have a widespread application to many types of rare transitional events, here only recent applications involving biomolecules are reviewed, including isomerization, protein folding, and enzyme catalysis.

The stochastic separatrix and the reaction coordinate for complex systems

We present a new approach to the identification of degrees of freedom which comprise a reaction coordinate in a complex system. The method begins with the generation of an ensemble of reactive trajectories. Each trajectory is analyzed for its equicommittor position or transition state; then the transition state ensemble is identified as the stochastic separatrix. Numerical analysis of the points along the separatrix for variability of coordinate location correctly identifies the components of the reaction coordinate in a test system of a double well coupled to a promoting vibration and a bath of linearly coupled oscillators.

Analysis of the free-energy surface of proteins from reversible folding simulations

Computer generated trajectories can, in principle, reveal the folding pathways of a protein at atomic resolution and possibly suggest general and simple rules for predicting the folded structure of a given sequence. While such reversible folding trajectories can only be determined ab initio using all-atom transferable force-fields for a few small proteins, they can be determined for a large number of proteins using coarse-grained and structure-based force-fields, in which a known folded structure is by construction the absolute energy and free-energy minimum. Here we use a model of the fast folding helical λ-repressor protein to generate trajectories in which native and non-native states are in equilibrium and transitions are accurately sampled. Yet, representation of the free-energy surface, which underlies the thermodynamic and dynamic properties of the protein model, from such a trajectory remains a challenge. Projections over one or a small number of arbitrarily chosen progress variables often hide the most important features of such surfaces. The results unequivocally show that an unprojected representation of the free-energy surface provides important and unbiased information and allows a simple and meaningful description of many-dimensional, heterogeneous trajectories, providing new insight into the possible mechanisms of fast-folding proteins.

The process of protein folding is a complex transition from a disordered to an ordered state. Here, we simulate a specific fast-folding protein at the point at which the native and denatured states are at equilibrium and show that obtaining an accurate description of the mechanisms of folding and unfolding is far from trivial. Using simple quantities which quantify the degree of native order is, in the case of this protein, clearly misleading. We show that an unbiased representation of the free-energy surface can be obtained; using such a representation we are able to redesign the landscape and thus modify, upon site-specific “mutations”, the folding and unfolding rates. This leads us to formulate a hypothesis to explain the very fast folding of many proteins.

Allosteric transitions of supramolecular systems explored by network models: application to chaperonin GroEL

Identification of pathways involved in the structural transitions of biomolecular systems is often complicated by the transient nature of the conformations visited across energy barriers and the multiplicity of paths accessible in the multidimensional energy landscape. This task becomes even more challenging in exploring molecular systems on the order of megadaltons. Coarse-grained models that lend themselves to analytical solutions appear to be the only possible means of approaching such cases. Motivated by the utility of elastic network models for describing the collective dynamics of biomolecular systems and by the growing theoretical and experimental evidence in support of the intrinsic accessibility of functional substates, we introduce a new method, adaptive anisotropic network model (aANM), for exploring functional transitions. Application to bacterial chaperonin GroEL and comparisons with experimental data, results from action minimization algorithm, and previous simulations support the utility of aANM as a computationally efficient, yet physically plausible, tool for unraveling potential transition pathways sampled by large complexes/assemblies. An important outcome is the assessment of the critical inter-residue interactions formed/broken near the transition state(s), most of which involve conserved residues.

Most proteins are biomolecular machines. They perform their function by undergoing changes between different structures. Understanding the mechanism of transition between these structures is of major importance to design methods for controlling such transitions, and thereby modulating protein function. Although there are many computational methods for exploring the transitions of small proteins, the task of exploring the transition pathways becomes prohibitively expensive in the case of supramolecular systems. The bacterial chaperonin GroEL is such a supramolecular machine. It plays an important role in assisting protein folding. During its function, GroEL undergoes structural transitions between multiple forms. Here, we are introducing a new methodology, based on elastic network models, for elucidating the transition mechanisms in such supramolecular systems. Application to GroEL provides us with biologically significant information on critical interactions and sequence of events occurring during the chaperonin machinery and key contacts that make and break at the transition. The method can be readily applied to other supramolecular machines.

Probing the nanosecond dynamics of a designed three-stranded beta-sheet with a massively parallel molecular dynamics simulation

Recently a temperature-jump FTIR study of a designed three-stranded sheet showing a fast relaxation time of approximately 140 +/- 20 ns was published. We performed massively parallel molecular dynamics simulations in explicit solvent to probe the structural events involved in this relaxation. While our simulations produce similar relaxation rates, the structural ensemble is broad. We observe the formation of turn structure, but only very weak interaction in the strand regions, which is consistent with the lack of strong backbone-backbone NOEs in previous structural NMR studies. These results suggest that either (D)P(D)P-II folds at time scales longer than 240 ns, or that (D)P(D)P-II is not a well-defined three-stranded beta-sheet. This work also provides an opportunity to compare the performance of several popular forcefield models against one another.

Thermodynamics and kinetics of protein folding under confinement

Understanding the effects of confinement on protein stability and folding kinetics is important for describing protein folding in the cellular environment. We have investigated the effects of confinement on two structurally distinct proteins as a function of the dimension d(c) and characteristic size R of the confining boundary. We find that the stabilization of the folded state relative to bulk conditions is quantitatively described by R(-gamma(c)), where the exponent gamma(c) is approximately 5/3 independent of the dimension of confinement d(c) (cylindrical, planar, or spherical). Moreover, we find that the logarithm of the folding rates also scale as R(-gamma(c)), with deviations only being seen for very small confining geometries, where folding is downhill; for both stability and kinetics, the dominant effect is the change in the free energy of the unfolded state. A secondary effect on the kinetics is a slight destabilization of the transition state by confinement, although the contacts present in the confined transition state are essentially identical to the bulk case. We investigate the effect of confinement on the position-dependent diffusion coefficients D(Q) for dynamics along the reaction coordinate Q (fraction of native contacts). The diffusion coefficients only change in the unfolded state basin, where they are increased because of compaction.

Chemical, physical, and theoretical kinetics of an ultrafast folding protein

An extensive set of equilibrium and kinetic data is presented and analyzed for an ultrafast folding protein--the villin subdomain. The equilibrium data consist of the excess heat capacity, tryptophan fluorescence quantum yield, and natural circular-dichroism spectrum as a function of temperature, and the kinetic data consist of time courses of the quantum yield from nanosecond-laser temperature-jump experiments. The data are well fit with three kinds of models--a three-state chemical-kinetics model, a physical-kinetics model, and an Ising-like theoretical model that considers 10(5) possible conformations (microstates). In both the physical-kinetics and theoretical models, folding is described as diffusion on a one-dimensional free-energy surface. In the physical-kinetics model the reaction coordinate is unspecified, whereas in the theoretical model, order parameters, either the fraction of native contacts or the number of native residues, are used as reaction coordinates. The validity of these two reaction coordinates is demonstrated from calculation of the splitting probability from the rate matrix of the master equation for all 10(5) microstates. The analysis of the data on site-directed mutants using the chemical-kinetics model provides information on the structure of the transition-state ensemble; the physical-kinetics model allows an estimate of the height of the free-energy barrier separating the folded and unfolded states; and the theoretical model provides a detailed picture of the free-energy surface and a residue-by-residue description of the evolution of the folded structure, yet contains many fewer adjustable parameters than either the chemical- or physical-kinetics models.

Theory, analysis, and interpretation of single-molecule force spectroscopy experiments

Dynamic force spectroscopy probes the kinetic and thermodynamic properties of single molecules and molecular assemblies. Here, we propose a simple procedure to extract kinetic information from such experiments. The cornerstone of our method is a transformation of the rupture-force histograms obtained at different force-loading rates into the force-dependent lifetimes measurable in constant-force experiments. To interpret the force-dependent lifetimes, we derive a generalization of Bell's formula that is formally exact within the framework of Kramers theory. This result complements the analytical expression for the lifetime that we derived previously for a class of model potentials. We illustrate our procedure by analyzing the nanopore unzipping of DNA hairpins and the unfolding of a protein attached by flexible linkers to an atomic force microscope. Our procedure to transform rupture-force histograms into the force-dependent lifetimes remains valid even when the molecular extension is a poor reaction coordinate and higher-dimensional free-energy surfaces must be considered. In this case the microscopic interpretation of the lifetimes becomes more challenging because the lifetimes can reveal richer, and even nonmonotonic, dependence on the force.

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.

Diffusive reaction dynamics on invariant free energy profiles

A fundamental problem in the analysis of protein folding and other complex reactions in which the entropy plays an important role is the determination of the activation free energy from experimental measurements or computer simulations. This article shows how to combine minimum-cut-based free energy profiles (F(C)), obtained from equilibrium molecular dynamics simulations, with conventional histogram-based free energy profiles (F(H)) to extract the coordinate-dependent diffusion coefficient on the F(C) (i.e., the method determines free energies and a diffusive preexponential factor along an appropriate reaction coordinate). The F(C), in contrast to the F(H), is shown to be invariant with respect to arbitrary transformations of the reaction coordinate, which makes possible partition of configuration space into basins in an invariant way. A "natural coordinate," for which F(H) and F(C) differ by a multiplicative constant (constant diffusion coefficient), is introduced. The approach is illustrated by a model one-dimensional system, the alanine dipeptide, and the folding reaction of a double beta-hairpin miniprotein. It is shown how the results can be used to test whether the putative reaction coordinate is a good reaction coordinate.

Macroscopically ordered water in nanopores

Water confined into the interior channels of narrow carbon nanotubes or transmembrane proteins forms collectively oriented molecular wires held together by tight hydrogen bonds. Here, we explore the thermodynamic stability and dipolar orientation of such 1D water chains from nanoscopic to macroscopic dimensions. We show that a dipole lattice model accurately recovers key properties of 1D confined water when compared to atomically detailed simulations. In a major reduction in computational complexity, we represent the dipole model in terms of effective Coulombic charges, which allows us to study pores of macroscopic lengths in equilibrium with a water bath (or vapor). We find that at ambient conditions, the water chains filling the tube are essentially continuous up to macroscopic dimensions. At reduced water vapor pressure, we observe a 1D Ising-like filling/emptying transition without a true phase transition in the thermodynamic limit. In the filled state, the chains of water molecules in the tube remain dipole-ordered up to macroscopic lengths of approximately 0.1 mm, and the dipolar order is estimated to persist for times up to approximately 0.1 s. The observed dipolar order in continuous water chains is a precondition for the use of nanoconfined 1D water as mediator of fast long-range proton transport, e.g., in fuel cells. For water-filled nanotube bundles and membranes, we expect anti-ferroelectric behavior, resulting in a rich phase diagram similar to that of a 2D Coulomb gas.

One-dimensional barrier-preserving free-energy projections of a beta-sheet miniprotein: new insights into the folding process

The conformational space of a 20-residue three-stranded antiparallel beta-sheet peptide (double hairpin) was sampled by equilibrium folding/unfolding molecular dynamics simulations for a total of 20 micros. The resulting one-dimensional free-energy profiles (FEPs) provide a detailed description of the free-energy basins and barriers for the folding reaction. The similarity of the FEPs obtained using the probability of folding before unfolding (pfold) or the mean first passage time supports the robustness of the procedure. The folded state and the most populated free-energy basins in the denatured state are described by the one-dimensional FEPs, which avoid the overlap of states present in the usual one- or two-dimensional projections. Within the denatured state, a basin with fluctuating helical conformations and a heterogeneous entropic state are populated near the melting temperature at about 11% and 33%, respectively. Folding pathways from the helical basin or enthalpic traps (with only one of the two hairpins formed) reach the native state through the entropic state, which is on-pathway and is separated by a low barrier from the folded state. A simplified equilibrium kinetic network based on the FEPs shows the complexity of the folding reaction and indicates, as augmented by additional analyses, that the basins in the denatured state are connected primarily by the native state. The overall folding kinetics shows single-exponential behavior because barriers between the non-native basins and the folded state have similar heights.

Normal mode partitioning of Langevin dynamics for biomolecules

We propose a novel normal mode multiple time stepping Langevin dynamics integrator called NML. The aim is to approximate the kinetics or thermodynamics of a biomolecule by a reduced model based on a normal mode decomposition of the dynamical space. Our basis set uses the eigenvectors of a mass reweighted Hessian matrix calculated with a biomolecular force field. This particular choice has the advantage of an ordering according to the eigenvalues, which have a physical meaning of being the square of the mode frequency. Low frequency eigenvalues correspond to more collective motions, whereas the highest frequency eigenvalues are the limiting factor for the stability of the integrator. In NML, the higher frequency modes are overdamped and relaxed near their energy minimum while respecting the subspace of low frequency dynamical modes. Our numerical results confirm that both sampling and rates are conserved for an implicitly solvated alanine dipeptide model, with only 30% of the modes propagated, when compared to the full model. For implicitly solvated systems, NML gives a twofold improvement in efficiency over plain Langevin dynamics for sampling a small 22 atom (alanine dipeptide) model and in excess of an order of magnitude for sampling an 882 atom (bovine pancreatic trypsin inhibitor) model, with good scaling with system size subject to the number of modes propagated. NML has been implemented in the open source software PROTOMOL.

Strongly anisotropic orientational relaxation of water molecules in narrow carbon nanotubes and nanorings

Molecular dynamics simulations of the orientational dynamics of water molecules confined in narrow carbon nanotubes and nanorings reveal that confinement leads to strong anisotropy in the orientational relaxation. The relaxation of the aligned dipole moments, occurring on a time scale of nanoseconds, is 3 orders of magnitude slower than that of bulk water. In contrast, the relaxation of the vector joining the two hydrogens is ten times faster compared to bulk, with a time scale of about 150 fs. The slow dipolar relaxation is mediated by the hopping of orientational defects, which are nucleated by the water molecules outside the tube, across the linear water chain.

Binding-induced folding of a natively unstructured transcription factor

Transcription factors are central components of the intracellular regulatory networks that control gene expression. An increasingly recognized phenomenon among human transcription factors is the formation of structure upon target binding. Here, we study the folding and binding of the pKID domain of CREB to the KIX domain of the co-activator CBP. Our simulations of a topology-based Gō-type model predict a coupled folding and binding mechanism, and the existence of partially bound intermediates. From transition-path and Phi-value analyses, we find that the binding transition state resembles the unstructured state in solution, implying that CREB becomes structured only after committing to binding. A change of structure following binding is reminiscent of an induced-fit mechanism and contrasts with models in which binding occurs to pre-structured conformations that exist in the unbound state at equilibrium. Interestingly, increasing the amount of structure in the unbound pKID reduces the rate of binding, suggesting a "fly-casting"-like process. We find that the inclusion of attractive non-native interactions results in the formation of non-specific encounter complexes that enhance the on-rate of binding, but do not significantly change the binding mechanism. Our study helps explain how being unstructured can confer an advantage in protein target recognition. The simulations are in general agreement with the results of a recently reported nuclear magnetic resonance study, and aid in the interpretation of the experimental binding kinetics.

Finding transition pathways using the string method with swarms of trajectories

An approach to find transition pathways in complex systems is presented. The method, which is related to the string method in collective variables of Maragliano et al. (J. Chem. Phys. 2006, 125, 024106), is conceptually simple and straightforward to implement. It consists of refining a putative transition path in the multidimensional space supported by a set of collective variables using the average dynamic drift of those variables. This drift is estimated on-the-fly via swarms of short unbiased trajectories started at different points along the path. Successive iterations of this algorithm, which can be naturally distributed over many computer nodes with negligible interprocessor communication, refine an initial trial path toward the most probable transition path (MPTP) between two stable basins. The method is first tested by determining the pathway for the C7eq to C7ax transition in an all-atom model of the alanine dipeptide in vacuum, which has been studied previously with the string method in collective variables. A transition path is found with a committor distribution peaked at 1/2 near the free energy maximum, in accord with previous results. Last, the method is applied to the allosteric conformational change in the nitrogen regulatory protein C (NtrC), represented here with a two-state elastic network model. Even though more than 550 collective variables are used to describe the conformational change, the path converges rapidly. Again, the committor distribution is found to be peaked around 1/2 near the free energy maximum between the two stable states, confirming that a genuine transition state has been localized in this complex multidimensional system.

Water permeation across nanochannels with defects

Defects are common in fabricated nanochannels. In this paper, water permeation across a single-walled carbon nanotube with defects was studied using molecular dynamics simulations. It is found that the impact on water permeation is negligible when the density of the defects is small, while a significant reduction in water permeation is observed when the density of the defects is high. These findings should be helpful in both understanding water permeation across nanochannels and designing efficient artificial nanochannel.

Solvent viscosity dependence of the protein folding dynamics

Solvent viscosity has been frequently adopted as an adjustable parameter in various computational studies (e.g., protein folding simulations) with implicit solvent models. A common approach is to use low viscosities to expedite simulations. While using viscosities lower than that of aqueous is unphysical, such treatment is based on observations that the viscosity affects the kinetics (rates) in a well-defined manner as described by Kramers' theory. Here, we investigate the effect of viscosity on the detailed dynamics (mechanism) of protein folding. On the basis of a simple mathematical model, we first show that viscosity may indeed affect the dynamics in a complex way. By applying the model to the folding of a small protein, we demonstrate that the detailed dynamics is affected rather pronouncedly especially at unphysically low viscosities, cautioning against using such viscosities. In this regard, our model may also serve as a diagnostic tool for validating low-viscosity simulations. It is also suggested that the viscosity dependence can be further exploited to gain information about the protein folding mechanism.

The protein folding problem

The "protein folding problem" consists of three closely related puzzles: (a) What is the folding code? (b) What is the folding mechanism? (c) Can we predict the native structure of a protein from its amino acid sequence? Once regarded as a grand challenge, protein folding has seen great progress in recent years. Now, foldable proteins and nonbiological polymers are being designed routinely and moving toward successful applications. The structures of small proteins are now often well predicted by computer methods. And, there is now a testable explanation for how a protein can fold so quickly: A protein solves its large global optimization problem as a series of smaller local optimization problems, growing and assembling the native structure from peptide fragments, local structures first.

Extensions to the likelihood maximization approach for finding reaction coordinates

This paper extends our previous work on obtaining reaction coordinates from aimless shooting and likelihood maximization. We introduce a simplified version of aimless shooting and a half-trajectory likelihood score based on the committor probability. Additionally, we analyze and compare the absolute log-likelihood score for perfect and approximate reaction coordinates. We also compare the aimless shooting and likelihood maximization approach to the earlier genetic neural network (GNN) approach of Ma and Dinner [J. Phys. Chem. B 109, 6769 (2005)]. For a fixed number of total trajectories in the GNN approach, the accuracy of the transition state ensemble decreases as the number of trajectories per committor probability estimate increases. This quantitatively demonstrates the benefit of individual committor probability realizations over committor probability estimates. Furthermore, when the least squares score of the GNN approach is applied to individual committor probability realizations, the likelihood score still provides a better approximation to the true transition state surface. Finally, the polymorph transition in terephthalic acid demonstrates that the new half-trajectory likelihood scheme estimates the transition state location more accurately than likelihood schemes based on the probability of being on a transition path.

Reorientional dynamics of water molecules in anionic hydration shells

Water molecule rotational dynamics within a chloride anion's first hydration shell are investigated through simulations. In contrast to recent suggestions that the ion's hydration shell is rigid during a water's reorientation, we find a labile hydration sphere, consistent with previous assessments of chloride as a weak structure breaker. The nondiffusive reorientation mechanism found involves a hydrogen-bond partner switch with a large amplitude angular jump and the water's departure from the anion's shell. An analytic extended jump model accounts for the simulation results, as well as available NMR and ultrafast spectroscopic data, and resolves the discrepancy between them.

Folding and unfolding of a photoswitchable peptide from picoseconds to microseconds

Using time-resolved IR spectroscopy, we monitored the kinetics of folding and unfolding processes of a photoswitchable 16-residue alanine-based alpha-helical peptide on a timescale from few picoseconds to almost 40 micros and over a large temperature range (279-318 K). The folding and unfolding processes were triggered by an ultrafast laser pulse that isomerized the cross linker within a few picoseconds. The main folding and unfolding times (700 ns and 150 ns, respectively, at room temperature) are in line with previous T-jump experiments obtained from similar peptides. However, both processes show complex, strongly temperature-dependent spectral kinetics that deviate clearly from a single-exponential behavior. Whereas in the unfolding experiment the ensemble starts from a well defined folded state, the starting ensemble in the folding experiment is more heterogeneous, which leads to distinctly different kinetics of the experiments, because they are sensitive to different regions of the energy surface. A qualitative agreement with the experimental data-set can be obtained by a model where the unfolded states act as a hub connected to several separated "misfolded" states with a distribution of rates. We conclude that a rather large spread of rates (k(1) : k(n) approximately 9) is needed to explain the experimentally observed stretched exponential response with stretching factor beta = 0.8 at 279 K.

Using the histogram test to quantify reaction coordinate error

Many schemes for calculating reaction rates and free energy barriers require an accurate reaction coordinate, but it is difficult to quantify reaction coordinate accuracy for complex processes like protein folding and nucleation. The histogram test, based on estimated committor probabilities, is often used as a qualitative indicator for good reaction coordinates. This paper derives the mean and variance of the intrinsic committor distribution in terms of the mean and variance of the histogram of committor estimates. These convenient formulas enable the first quantitative calculations of reaction coordinate error for complex systems. An example shows that the approximate transition state surface from Peters' and Trout's reaction coordinate for nucleation in the Ising model gives a mean committor probability of 0.495 and a standard deviation of 0.042.