Finite Temperature String Method for the Study of Rare Events†

Applications of rare event dynamics on the free energy calculations for membrane protein systems

Techniques of rare event dynamics were reviewed including the string methods, which will be implemented with the biochemical simulation packages. The existing methods were applied to study biological systems with relevance to drug design and drug metabolism. The rare event dynamics simulations were performed to understand the kinetic and thermodynamic free energy information on the drug binding sites in the M2 proton channel, and also the free energy of insertion and association of membrane proteins and membrane active peptides. Results give a theoretical framework to interpret and reconcile existing and often conflicting opinions.

Efficiently computing pathway free energies: New approaches based on chain-of-replica and Non-Boltzmann Bennett reweighting schemes

BACKGROUND

Accurately modeling condensed phase processes is one of computation's most difficult challenges. Include the possibility that conformational dynamics may be coupled to chemical reactions, where multiscale (i.e., QM/MM) methods are needed, and this task becomes even more daunting.

METHODS

Free energy simulations (i.e., molecular dynamics), multiscale modeling, and reweighting schemes.

RESULTS

Herein, we present two new approaches for mitigating the aforementioned challenges. The first is a new chain-of-replica method (off-path simulations, OPS) for computing potentials of mean force (PMFs) along an easily defined reaction coordinate. This development is coupled with a new distributed, highly-parallel replica framework (REPDstr) within the CHARMM package. Validation of these new schemes is carried out on two processes that undergo conformational changes. First is the simple torsional rotation of butane, while a much more challenging glycosidic rotation (in vacuo and solvated) is the second. Additionally, a new approach that greatly improves (i.e., possibly an order of magnitude) the efficiency of computing QM/MM PMFs is introduced and compared to standard schemes. Our efforts are grounded in the recently developed method for efficiently computing QM-based free energies (i.e., QM-Non-Boltzmann Bennett, QM-NBB). Again, we validate this new technique by computing the QM/MM PMF of butane's torsional rotation.

CONCLUSIONS

The OPS-REPDstr method is a promising new approach that overcomes many limitations of standard pathway simulations in CHARMM. The combination of QM-NBB with pathway techniques is very promising as it offers significant advantages over current procedures.

GENERAL SIGNIFICANCE

Efficiently computing potentials of mean force is a major, unresolved, area of interest. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

Bézier curve string method for the study of rare events in complex chemical systems

We present a new string method for finding the most probable transition pathway and optimal reaction coordinate in complex chemical systems. Our approach evolves an analytic parametric curve, known as a Bézier curve, to the most probable transition path between metastable regions in configuration space. In addition, we demonstrate that the geometric properties of the Bézier curve can be used to construct the optimal reaction coordinate near the most probable reaction path, and can further be used to devise a ranking vector capable of identifying precisely which collective variables are most important for governing the transition between metastable states. We discuss the algorithmic details of the Bézier curve string method, analyze its stability, accuracy and efficiency, and illustrate its capabilities using model potential energy functions. In particular, we use the degree elevation property of Bézier curves to develop an algorithm that adaptively learns the degree polynomial necessary to accurately represent the most probable transition path. Subsequently, we apply our method to the isomerization of alanine dipeptide, and demonstrate that the reaction coordinate obtained from the Bézier curve string method is in excellent agreement with the optimal reaction coordinate constructed from an aimless shooting and maximum likelihood procedure. Finally, we apply our method to a large complex system and study the homogenous nucleation of benzene from the melt. In these two examples, we illustrate that the ranking vector correctly identifies which collective variables govern these chemical transitions.

Towards fast, rigorous and efficient conformational sampling of biomolecules: Advances in accelerated molecular dynamics

BACKGROUND

Accelerated molecular dynamics (aMD) has been proven to be a powerful biasing method for enhanced sampling of biomolecular conformations on general-purpose computational platforms. Biologically important long timescale events that are beyond the reach of standard molecular dynamics can be accessed without losing the detailed atomistic description of the system in aMD. Over other biasing methods, aMD offers the advantages of tuning the level of acceleration to access the desired timescale without any advance knowledge of the reaction coordinate.

SCOPE OF REVIEW

Recent advances in the implementation of aMD and its applications to small peptides and biological macromolecules are reviewed here along with a brief account of all the aMD variants introduced in the last decade.

MAJOR CONCLUSIONS

In comparison to the original implementation of aMD, the recent variant in which all the rotatable dihedral angles are accelerated (RaMD) exhibits faster convergence rates and significant improvement in statistical accuracy of retrieved thermodynamic properties. RaMD in conjunction with accelerating diffusive degrees of freedom, i.e. dual boosting, has been rigorously tested for the most difficult conformational sampling problem, protein folding. It has been shown that RaMD with dual boosting is capable of efficiently sampling multiple folding and unfolding events in small fast folding proteins.

GENERAL SIGNIFICANCE

RaMD with the dual boost approach opens exciting possibilities for sampling multiple timescales in biomolecules. While equilibrium properties can be recovered satisfactorily from aMD-based methods, directly obtaining dynamics and kinetic rates for larger systems presents a future challenge. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

Modeling and enhanced sampling of molecular systems with smooth and nonlinear data-driven collective variables

Collective variables (CVs) are low-dimensional representations of the state of a complex system, which help us rationalize molecular conformations and sample free energy landscapes with molecular dynamics simulations. Given their importance, there is need for systematic methods that effectively identify CVs for complex systems. In recent years, nonlinear manifold learning has shown its ability to automatically characterize molecular collective behavior. Unfortunately, these methods fail to provide a differentiable function mapping high-dimensional configurations to their low-dimensional representation, as required in enhanced sampling methods. We introduce a methodology that, starting from an ensemble representative of molecular flexibility, builds smooth and nonlinear data-driven collective variables (SandCV) from the output of nonlinear manifold learning algorithms. We demonstrate the method with a standard benchmark molecule, alanine dipeptide, and show how it can be non-intrusively combined with off-the-shelf enhanced sampling methods, here the adaptive biasing force method. We illustrate how enhanced sampling simulations with SandCV can explore regions that were poorly sampled in the original molecular ensemble. We further explore the transferability of SandCV from a simpler system, alanine dipeptide in vacuum, to a more complex system, alanine dipeptide in explicit water.

Reaction Path Force Matching: A New Strategy of Fitting Specific Reaction Parameters for Semiempirical Methods in Combined QM/MM Simulations

We present a general strategy of reparametrizing semiempirical (SE) methods against ab initio (AI) methods for combined quantum mechanical and molecular mechanical (QM/MM) simulations of specific chemical reactions in condensed phases. The resulting approach, designated Reaction Path Force Matching (RP-FM), features cycles of sampling configurations along a reaction path on an efficient SE/MM potential energy surface (PES) and adjusting specific reaction parameters (SRPs) in the SE method such that the atomic forces computed at the target AI/MM level are reproduced. Iterative applications of the RP-FM cycle make possible achieving the accuracy of AI/MM simulations without explicitly sampling the computationally expensive AI/MM PES. The bypassed sampling, nevertheless, is implicitly accomplished through the aid of the efficient SE-SRP/MM PES, on which the target-level reaction path is expected to be obtained upon convergence. We demonstrate the effectiveness of the RP-FM procedure for a symmetric proton transfer reaction in the gas phase and in solution. The remarkable agreements between the RP-FM optimized SE-SRP methods and the target AI method on various properties, including energy profiles, potential of mean force free energy profiles, atomic forces, charge populations, and solvation effects, suggest that RP-FM can be used as an efficient and reliable strategy for simulating condensed-phase chemical reactions.

Using multiscale preconditioning to accelerate the convergence of iterative molecular calculations

Iterative procedures for optimizing properties of molecular models often converge slowly owing to the computational cost of accurately representing features of interest. Here, we introduce a preconditioning scheme that allows one to use a less expensive model to guide exploration of the energy landscape of a more expensive model and thus speed the discovery of locally stable states of the latter. We illustrate our approach in the contexts of energy minimization and the string method for finding transition pathways. The relation of the method to other multilevel simulation techniques and possible extensions are discussed.

Molecular recognition of DNA by ligands: roughness and complexity of the free energy profile

Understanding the molecular mechanism by which probes and chemotherapeutic agents bind to nucleic acids is a fundamental issue in modern drug design. From a computational perspective, valuable insights are gained by the estimation of free energy landscapes as a function of some collective variables (CVs), which are associated with the molecular recognition event. Unfortunately the choice of CVs is highly non-trivial because of DNA's high flexibility and the presence of multiple association-dissociation events at different locations and/or sliding within the grooves. Here we have applied a modified version of Locally-Scaled Diffusion Map (LSDMap), a nonlinear dimensionality reduction technique for decoupling multiple-timescale dynamics in macromolecular systems, to a metadynamics-based free energy landscape calculated using a set of intuitive CVs. We investigated the binding of the organic drug anthramycin to a DNA 14-mer duplex. By performing an extensive set of metadynamics simulations, we observed sliding of anthramycin along the full-length DNA minor groove, as well as several detachments from multiple sites, including the one identified by X-ray crystallography. As in the case of equilibrium processes, the LSDMap analysis is able to extract the most relevant collective motions, which are associated with the slow processes within the system, i.e., ligand diffusion along the minor groove and dissociation from it. Thus, LSDMap in combination with metadynamics (and possibly every equivalent method) emerges as a powerful method to describe the energetics of ligand binding to DNA without resorting to intuitive ad hoc reaction coordinates.

Generalized Scalable Multiple Copy Algorithms for Molecular Dynamics Simulations in NAMD

Computational methodologies that couple the dynamical evolution of a set of replicated copies of a system of interest offer powerful and flexible approaches to characterize complex molecular processes. Such multiple copy algorithms (MCAs) can be used to enhance sampling, compute reversible work and free energies, as well as refine transition pathways. Widely used examples of MCAs include temperature and Hamiltonian-tempering replica-exchange molecular dynamics (T-REMD and H-REMD), alchemical free energy perturbation with lambda replica-exchange (FEP/λ-REMD), umbrella sampling with Hamiltonian replica exchange (US/H-REMD), and string method with swarms-of-trajectories conformational transition pathways. Here, we report a robust and general implementation of MCAs for molecular dynamics (MD) simulations in the highly scalable program NAMD built upon the parallel programming system Charm++. Multiple concurrent NAMD instances are launched with internal partitions of Charm++ and located continuously within a single communication world. Messages between NAMD instances are passed by low-level point-to-point communication functions, which are accessible through NAMD's Tcl scripting interface. The communication-enabled Tcl scripting provides a sustainable application interface for end users to realize generalized MCAs without modifying the source code. Illustrative applications of MCAs with fine-grained inter-copy communication structure, including global lambda exchange in FEP/λ-REMD, window swapping US/H-REMD in multidimensional order parameter space, and string method with swarms-of-trajectories were carried out on IBM Blue Gene/Q to demonstrate the versatility and massive scalability of the present implementation.

Calcium inhibition of ribonuclease H1 two-metal ion catalysis

Most phosphate-processing enzymes require Mg(2+) as a cofactor to catalyze nucleotide cleavage and transfer reactions. Ca(2+) ions inhibit many of these enzymatic activities, despite Ca(2+) and Mg(2+) having comparable binding affinities and overall biological abundances. Here we study the molecular details of the calcium inhibition mechanism for phosphodiester cleavage, an essential reaction in the metabolism of nucleic acids and nucleotides, by comparing Ca(2+)- and Mg(2+) catalyzed reactions. We study the functional roles of the specific metal ion sites A and B in enabling the catalytic cleavage of an RNA/DNA hybrid substrate by B. halodurans ribonuclease (RNase) H1 using hybrid quantum-mechanics/molecular mechanics (QM/MM) free energy calculations. We find that Ca(2+) substitution of either of the two active-site Mg(2+) ions substantially increases the height of the reaction barrier and thereby abolishes the catalytic activity. Remarkably, Ca(2+) at the A site is inactive also in Mg(2+)-optimized active-site structures along the reaction path, whereas Mg(2+) substitution recovers activity in Ca(2+)-optimized structures. Geometric changes resulting from Ca(2+) substitution at metal ion site A may thus be a secondary factor in the loss of catalytic activity. By contrast, at metal ion site B geometry plays a more important role, with only a partial recovery of activity after Mg(2+) substitution in Ca(2+)-optimized structures. Ca(2+)-substitution also leads to a change in mechanism, with deprotonation of the water nucleophile requiring a closer approach to the scissile phosphate, which in turn increases the barrier. As a result, Ca(2+) is less efficient in activating the water. As a likely cause for the different reactivities of Mg(2+) and Ca(2+) ions in site A, we identify differences in charge transfer to the ions and the associated decrease in the pKa of the oxygen nucleophile attacking the phosphate group.

Unbinding free energy of acetylcholinesterase bound oxime drugs along the gorge pathway from metadynamics-umbrella sampling investigation

Because of the pivotal role that the nerve enzyme, acetylcholinesterase plays in terminating nerve impulses at cholinergic synapses. Its active site, located deep inside a 20 Å gorge, is a vulnerable target of the lethal organophosphorus compounds. Potent reactivators of the intoxicated enzyme are nucleophiles, such as bispyridinium oxime that binds to the peripheral anionic site and the active site of the enzyme through suitable cation-π interactions. Atomic scale molecular dynamics and free energy calculations in explicit water are used to study unbinding pathways of two oxime drugs (Ortho-7 and Obidoxime) from the gorge of the enzyme. The role of enzyme-drug cation-π interactions are explored with the metadynamics simulation. The metadynamics discovered potential of mean force (PMF) of the unbinding events is refined by the umbrella sampling (US) corrections. The bidimensional free energy landscape of the metadynamics runs are further subjected to finite temperature string analysis to obtain the transition tube connecting the minima and bottlenecks of the unbinding pathway. The PMF is also obtained from US simulations using the biasing potential constructed from the transition tube and are found to be consistent with the metadynamics-US corrected results. Although experimental structural data clearly shows analogous coordination of the two drugs inside the gorge in the bound state, the PMF of the drug trafficking along the gorge pathway point, within an equilibrium free energy context, to a multistep process that differs from one another. Routes, milestones and subtlety toward the unbinding pathway of the two oximes at finite temperature are identified.

WExplore: hierarchical exploration of high-dimensional spaces using the weighted ensemble algorithm

As most relevant motions in biomolecular systems are inaccessible to conventional molecular dynamics simulations, algorithms that enhance sampling of rare events are indispensable. Increasing interest in intrinsically disordered systems and the desire to target ensembles of protein conformations (rather than single structures) in drug development motivate the need for enhanced sampling algorithms that are not limited to "two-basin" problems, and can efficiently determine structural ensembles. For systems that are not well-studied, this must often be done with little or no information about the dynamics of interest. Here we present a novel strategy to determine structural ensembles that uses dynamically defined sampling regions that are organized in a hierarchical framework. It is based on the weighted ensemble algorithm, where an ensemble of copies of the system ("replicas") is directed to new regions of configuration space through merging and cloning operations. The sampling hierarchy allows for a large number of regions to be defined, while using only a small number of replicas that can be balanced over multiple length scales. We demonstrate this algorithm on two model systems that are analytically solvable and examine the 10-residue peptide chignolin in explicit solvent. The latter system is analyzed using a configuration space network, and novel hydrogen bonds are found that facilitate folding.

Comparison between Mean Forces and Swarms-of-Trajectories String Methods

The original formulation of the string method in collective variable space is compared with a recent variant called string method with swarms-of-trajectories. The assumptions made in the original method are revisited and the significance of the minimum free energy path (MFEP) is discussed in the context of reactive events. These assumptions are compared to those made in the string method with swarms-of-trajectories, and shown to be equivalent in a certain regime: in particular an expression for the path identified by the swarms-of-trajectories method is given and shown to be closely related to the MFEP. Finally, the algorithmic aspects of both methods are compared.

Quantum mechanical/molecular mechanical free energy simulations of the self-cleavage reaction in the hepatitis delta virus ribozyme

The hepatitis delta virus (HDV) ribozyme catalyzes a self-cleavage reaction using a combination of nucleobase and metal ion catalysis. Both divalent and monovalent ions can catalyze this reaction, although the rate is slower with monovalent ions alone. Herein, we use quantum mechanical/molecular mechanical (QM/MM) free energy simulations to investigate the mechanism of this ribozyme and to elucidate the roles of the catalytic metal ion. With Mg²⁺ at the catalytic site, the self-cleavage mechanism is observed to be concerted with a phosphorane-like transition state and a free energy barrier of ∼13 kcal/mol, consistent with free energy barrier values extrapolated from experimental studies. With Na⁺ at the catalytic site, the mechanism is observed to be sequential, passing through a phosphorane intermediate, with free energy barriers of 2–4 kcal/mol for both steps; moreover, proton transfer from the exocyclic amine of protonated C75 to the nonbridging oxygen of the scissile phosphate occurs to stabilize the phosphorane intermediate in the sequential mechanism. To explain the slower rate observed experimentally with monovalent ions, we hypothesize that the activation of the O2′ nucleophile by deprotonation and orientation is less favorable with Na⁺ ions than with Mg²⁺ ions. To explore this hypothesis, we experimentally measure the pK_a of O2′ by kinetic and NMR methods and find it to be lower in the presence of divalent ions rather than only monovalent ions. The combined theoretical and experimental results indicate that the catalytic Mg²⁺ ion may play three key roles: assisting in the activation of the O2′ nucleophile, acidifying the general acid C75, and stabilizing the nonbridging oxygen to prevent proton transfer to it.

On the study of crystal growth via interfacial analysis and string optimization

Mathematical ‘strings’ can be used with computer simulations and statistical mechanics to calculate the fraction of growth units and activation energies of flexible molecules present at the crystal–solution interface.

Computational study of peptide permeation through membrane: Searching for hidden slow variables

Atomically detailed molecular dynamics trajectories in conjunction with Milestoning are used to analyze the different contributions of coarse variables to the permeation process of a small peptide (N-acetyl-L-tryptophanamide, NATA) through a 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) membrane. The peptide reverses its overall orientation as it permeates through the biological bilayer. The large change in orientation is investigated explicitly but is shown to impact the free energy landscape and permeation time only moderately. Nevertheless, a significant difference in permeation properties of the two halves of the membrane suggests the presence of other hidden slow variables. We speculate, based on calculation of the potential of mean force, that a conformational transition of NATA makes significant contribution to these differences. Other candidates for hidden slow variables may include water permeation and collective motions of phospholipids.

Elucidating the ensemble of functionally-relevant transitions in protein systems with a robotics-inspired method

Background

Many proteins tune their biological function by transitioning between different functional states, effectively acting as dynamic molecular machines. Detailed structural characterization of transition trajectories is central to understanding the relationship between protein dynamics and function. Computational approaches that build on the Molecular Dynamics framework are in principle able to model transition trajectories at great detail but also at considerable computational cost. Methods that delay consideration of dynamics and focus instead on elucidating energetically-credible conformational paths connecting two functionally-relevant structures provide a complementary approach. Effective sampling-based path planning methods originating in robotics have been recently proposed to produce conformational paths. These methods largely model short peptides or address large proteins by simplifying conformational space.

Methods

We propose a robotics-inspired method that connects two given structures of a protein by sampling conformational paths. The method focuses on small- to medium-size proteins, efficiently modeling structural deformations through the use of the molecular fragment replacement technique. In particular, the method grows a tree in conformational space rooted at the start structure, steering the tree to a goal region defined around the goal structure. We investigate various bias schemes over a progress coordinate for balance between coverage of conformational space and progress towards the goal. A geometric projection layer promotes path diversity. A reactive temperature scheme allows sampling of rare paths that cross energy barriers.

Results and conclusions

Experiments are conducted on small- to medium-size proteins of length up to 214 amino acids and with multiple known functionally-relevant states, some of which are more than 13Å apart of each-other. Analysis reveals that the method effectively obtains conformational paths connecting structural states that are significantly different. A detailed analysis on the depth and breadth of the tree suggests that a soft global bias over the progress coordinate enhances sampling and results in higher path diversity. The explicit geometric projection layer that biases the exploration away from over-sampled regions further increases coverage, often improving proximity to the goal by forcing the exploration to find new paths. The reactive temperature scheme is shown effective in increasing path diversity, particularly in difficult structural transitions with known high-energy barriers.

Analyzing milestoning networks for molecular kinetics: definitions, algorithms, and examples

Network representations are becoming increasingly popular for analyzing kinetic data from techniques like Milestoning, Markov State Models, and Transition Path Theory. Mapping continuous phase space trajectories into a relatively small number of discrete states helps in visualization of the data and in dissecting complex dynamics to concrete mechanisms. However, not only are molecular networks derived from molecular dynamics simulations growing in number, they are also getting increasingly complex, owing partly to the growth in computer power that allows us to generate longer and better converged trajectories. The increased complexity of the networks makes simple interpretation and qualitative insight of the molecular systems more difficult to achieve. In this paper, we focus on various network representations of kinetic data and algorithms to identify important edges and pathways in these networks. The kinetic data can be local and partial (such as the value of rate coefficients between states) or an exact solution to kinetic equations for the entire system (such as the stationary flux between vertices). In particular, we focus on the Milestoning method that provides fluxes as the main output. We proposed Global Maximum Weight Pathways as a useful tool for analyzing molecular mechanism in Milestoning networks. A closely related definition was made in the context of Transition Path Theory. We consider three algorithms to find Global Maximum Weight Pathways: Recursive Dijkstra's, Edge-Elimination, and Edge-List Bisection. The asymptotic efficiency of the algorithms is analyzed and numerical tests on finite networks show that Edge-List Bisection and Recursive Dijkstra's algorithms are most efficient for sparse and dense networks, respectively. Pathways are illustrated for two examples: helix unfolding and membrane permeation. Finally, we illustrate that networks based on local kinetic information can lead to incorrect interpretation of molecular mechanisms.

Intermediate state trapping of a voltage sensor

Voltage sensor domains (VSDs) regulate ion channels and enzymes by undergoing conformational changes depending on membrane electrical signals. The molecular mechanisms underlying the VSD transitions are not fully understood. Here, we show that some mutations of I241 in the S1 segment of the Shaker Kv channel positively shift the voltage dependence of the VSD movement and alter the functional coupling between VSD and pore domains. Among the I241 mutants, I241W immobilized the VSD movement during activation and deactivation, approximately halfway between the resting and active states, and drastically shifted the voltage activation of the ionic conductance. This phenotype, which is consistent with a stabilization of an intermediate VSD conformation by the I241W mutation, was diminished by the charge-conserving R2K mutation but not by the charge-neutralizing R2Q mutation. Interestingly, most of these effects were reproduced by the F244W mutation located one helical turn above I241. Electrophysiology recordings using nonnatural indole derivatives ruled out the involvement of cation-Π interactions for the effects of the Trp inserted at positions I241 and F244 on the channel’s conductance, but showed that the indole nitrogen was important for the I241W phenotype. Insight into the molecular mechanisms responsible for the stabilization of the intermediate state were investigated by creating in silico the mutations I241W, I241W/R2K, and F244W in intermediate conformations obtained from a computational VSD transition pathway determined using the string method. The experimental results and computational analysis suggest that the phenotype of I241W may originate in the formation of a hydrogen bond between the indole nitrogen atom and the backbone carbonyl of R2. This work provides new information on intermediate states in voltage-gated ion channels with an approach that produces minimum chemical perturbation.

Insight into the molecular mechanism of water evaporation via the finite temperature string method

The process of water's evaporation at its liquid/air interface has proven challenging to study experimentally and, because it constitutes a rare event on molecular time scales, presents a challenge for computer simulations as well. In this work, we simulated water's evaporation using the classical extended simple point charge model water model, and identified a minimum free energy path for this process in terms of 10 descriptive order parameters. The measured free energy change was 7.4 kcal/mol at 298 K, in reasonable agreement with the experimental value of 6.3 kcal/mol, and the mean first-passage time was 1375 ns for a single molecule, corresponding to an evaporation coefficient of 0.25. In the observed minimum free energy process, the water molecule diffuses to the surface, and tends to rotate so that its dipole and one O-H bond are oriented outward as it crosses the Gibbs dividing surface. As the water molecule moves further outward through the interfacial region, its local density is higher than the time-averaged density, indicating a local solvation shell that protrudes from the interface. The water molecule loses donor and acceptor hydrogen bonds, and then, with its dipole nearly normal to the interface, stops donating its remaining hydrogen bond. At that point, when the final, accepted hydrogen bond is broken, the water molecule is free. We also analyzed which order parameters are most important in the process and in reactive trajectories, and found that the relative orientation of water molecules near the evaporating molecule, and the number of accepted hydrogen bonds, were important variables in reactive trajectories and in kinetic descriptions of the process.

Applications of rare event dynamics on the free energy calculations for membrane protein systems

Techniques of rare event dynamics were reviewed, including string methods, which will be implemented with the biochemical simulation packages. The existing methods were applied to study biological systems with relevance to drug design and drug metabolism. The rare event dynamics simulations were performed to understand the kinetic and thermodynamic free energy information on the drug binding sites in the M2 proton channel, the free energy of insertion and association of membrane proteins and membrane active peptides. Results give a theoretical framework to interpret and reconcile existing and often conflicting opinions.

Hidden Conformation Events in DNA Base Extrusions: A Generalized Ensemble Path Optimization and Equilibrium Simulation Study

DNA base extrusion is a crucial component of many biomolecular processes. Elucidating how bases are selectively extruded from the interiors of double-strand DNAs is pivotal to accurately understanding and efficiently sampling this general type of conformational transitions. In this work, the on-the-path random walk (OTPRW) method, which is the first generalized ensemble sampling scheme designed for finite-temperature-string path optimizations, was improved and applied to obtain the minimum free energy path (MFEP) and the free energy profile of a classical B-DNA major-groove base extrusion pathway. Along the MFEP, an intermediate state and the corresponding transition state were located and characterized. The MFEP result suggests that a base-plane-elongation event rather than the commonly focused base-flipping event is dominant in the transition state formation portion of the pathway; and the energetic penalty at the transition state is mainly introduced by the stretching of the Watson-Crick base pair. Moreover to facilitate the essential base-plane-elongation dynamics, the surrounding environment of the flipped base needs to be intimately involved. Further taking the advantage of the extended-dynamics nature of the OTPRW Hamiltonian, an equilibrium generalized ensemble simulation was performed along the optimized path; and based on the collected samples, several base-flipping (opening) angle collective variables were evaluated. In consistence with the MFEP result, the collective variable analysis result reveals that none of these commonly employed flipping (opening) angles alone can adequately represent the base extrusion pathway, especially in the pre-transition-state portion. As further revealed by the collective variable analysis, the base-pairing partner of the extrusion target undergoes a series of in-plane rotations to facilitate the base-plane-elongation dynamics. A base-plane rotation angle is identified to be a possible reaction coordinate to represent these in-plane rotations. Notably, these in-plane rotation motions may play a pivotal role in determining the base extrusion selectivity.

Conformational dynamics of a ligand-free adenylate kinase

Adenylate kinase (AdK) is a phosphoryl-transfer enzyme with important physiological functions. Based on a ligand-free open structure and a ligand-bound closed structure solved by crystallography, here we use molecular dynamics simulations to examine the stability and dynamics of AdK conformations in the absence of ligands. We first perform multiple simulations starting from the open or the closed structure, and observe their free evolutions during a simulation time of 100 or 200 nanoseconds. In all seven simulations starting from the open structure, AdK remained stable near the initial conformation. The eight simulations initiated from the closed structure, in contrast, exhibited large variation in the subsequent evolutions, with most (seven) undergoing large-scale spontaneous conformational changes and approaching or reaching the open state. To characterize the thermodynamics of the transition, we propose and apply a new sampling method that employs a series of restrained simulations to calculate a one-dimensional free energy along a curved pathway in the high-dimensional conformational space. Our calculated free energy profile features a single minimum at the open conformation, and indicates that the closed state, with a high (∼13 kcal/mol) free energy, is not metastable, consistent with the observed behaviors of the unrestrained simulations. Collectively, our simulations suggest that it is energetically unfavorable for the ligand-free AdK to access the closed conformation, and imply that ligand binding may precede the closure of the enzyme.

All-atom structural models of insulin binding to the insulin receptor in the presence of a tandem hormone-binding element

Insulin regulates blood glucose levels in higher organisms by binding to and activating insulin receptor (IR), a constitutively homodimeric glycoprotein of the receptor tyrosine kinase (RTK) superfamily. Therapeutic efforts in treating diabetes have been significantly impeded by the absence of structural information on the activated form of the insulin/IR complex. Mutagenesis and photo-crosslinking experiments and structural information on insulin and apo-IR strongly suggest that the dual-chain insulin molecule, unlike the related single-chain insulin-like growth factors, binds to IR in a very different conformation than what is displayed in storage forms of the hormone. In particular, hydrophobic residues buried in the core of the folded insulin molecule engage the receptor. There is also the possibility of plasticity in the receptor structure based on these data, which may in part be due to rearrangement of the so-called CT-peptide, a tandem hormone-binding element of IR. These possibilities provide opportunity for large-scale molecular modeling to contribute to our understanding of this system. Using various atomistic simulation approaches, we have constructed all-atom structural models of hormone/receptor complexes in the presence of CT in its crystallographic position and a thermodynamically favorable displaced position. In the "displaced-CT" complex, many more insulin-receptor contacts suggested by experiments are satisfied, and our simulations also suggest that R-insulin potentially represents the receptor-bound form of hormone. The results presented in this work have further implications for the design of receptor-specific agonists/antagonists.

Simulating rare events using a weighted ensemble-based string method

We introduce an extension to the weighted ensemble (WE) path sampling method to restrict sampling to a one-dimensional path through a high dimensional phase space. Our method, which is based on the finite-temperature string method, permits efficient sampling of both equilibrium and non-equilibrium systems. Sampling obtained from the WE method guides the adaptive refinement of a Voronoi tessellation of order parameter space, whose generating points, upon convergence, coincide with the principle reaction pathway. We demonstrate the application of this method to several simple, two-dimensional models of driven Brownian motion and to the conformational change of the nitrogen regulatory protein C receiver domain using an elastic network model. The simplicity of the two-dimensional models allows us to directly compare the efficiency of the WE method to conventional brute force simulations and other path sampling algorithms, while the example of protein conformational change demonstrates how the method can be used to efficiently study transitions in the space of many collective variables.

Molecular machines

Molecular machines (MM) are essential components of living cells. They conduct mechanical work, transport materials into and out of cells, assist in processing enzymatic reactions, and more. Their operations are frequently combined with significant conformational transitions. Computational studies of these conformational transitions and their coupling to molecular functions are discussed. It is argued that coarse descriptions of these molecules which are based on mass density and shape provide useful information on directions of action. It is further argued that MM are likely to have well focused and narrow reaction pathways. The proposal for such pathways is supported by evolutionary analyses of homologous machines. Finally, these observations are used to build atomically detailed models of these systems that are making the link from structure to functions (kinetics and thermodynamics). For that purpose enhanced sampling techniques are required.

Theoretical study of the mechanism of the hydride transfer between ferredoxin-NADP+ reductase and NADP+: the role of Tyr303

During photosynthesis, ferredoxin-NADP(+) reductase (FNR) catalyzes the electron transfer from ferredoxin to NADP(+) via its FAD cofactor. The final hydride transfer event between FNR and the nucleotide is a reversible process. Two different transient charge-transfer complexes form prior to and upon hydride transfer, FNR(rd)-NADP(+) and FNR(ox)-NADPH, regardless of the hydride transfer direction. Experimental structures of the FNR(ox):NADP(+) interaction have suggested a series of conformational rearrangements that might contribute to attaining the catalytically competent complex, but to date, no direct experimental information about the structure of this complex is available. Recently, a molecular dynamics (MD) theoretical approach was used to provide a putative organization of the active site that might represent a structure close to the transient catalytically competent interaction of Anabaena FNR with its coenzyme, NADP(+). Using this structure, we performed fully microscopic simulations of the hydride transfer processes between Anabaena FNR(rd)/FNR(ox) and NADP(+)/H, accounting also for the solvation. A dual-level QM/MM hybrid approach was used to describe the potential energy surface of the whole system. MD calculations using the finite-temperature string method combined with the WHAM method provided the potential of mean force for the hydride transfer processes. The results confirmed that the structural model of the reactants evolves to a catalytically competent transition state through very similar free energy barriers for both the forward and reverse reactions, in good agreement with the experimental hydride transfer rate constants reported for this system. This theoretical approach additionally provides subtle structural details of the mechanism in wild-type FNR and provides an explanation why Tyr303 makes possible the photosynthetic reaction, a process that cannot occur when this Tyr is replaced by a Ser.

Automated Transition State Searches without Evaluating the Hessian

Accurate and speedy determination of transition structures (TSs) is essential for computational studies on reaction pathways, particularly when the process involves expensive electronic structure calculations. Many search algorithms require a good initial guess of the TS geometry, as well as a Hessian input that possesses a structure consistent with the desired saddle point. Among the double-ended interpolation methods for generation of the guess for the TS, the freezing string method (FSM) is proven to be far less expensive compared to its predecessor, the growing string method (GSM). In this paper, it is demonstrated that the efficiency of this technique can be improved further by replacing the conjugate gradient optimization step (FSM-CG) with a quasi-Newton line search coupled with a BFGS Hessian update (FSM-BFGS). A second crucial factor that affects the speed with which convergence to the TS is achieved is the quality and cost of the Hessian of the energy for the guessed TS. For electronic structure calculations, the cost of calculating an exact Hessian increases more rapidly with system size than the energy and gradient. Therefore, to sidestep calculation of the exact Hessian, an approximate Hessian is constructed, using the tangent direction and local curvature at the TS guess. It is demonstrated that the partitioned-rational function optimization algorithm for locating TSs with this approximate Hessian input performs at least as well as with an exact Hessian input in most test cases. The two techniques, FSM and approximate Hessian construction, therefore can significantly reduce costs associated with finding TSs.

"DFG-flip" in the insulin receptor kinase is facilitated by a helical intermediate state of the activation loop

We have characterized a large-scale inactive-to-active conformational change in the activation-loop of the insulin receptor kinase domain at the atomistic level via untargeted temperature-accelerated molecular dynamics (TAMD) and free-energy calculations using the string method. TAMD simulations consistently show folding of the A-loop into a helical conformation followed by unfolding to an active conformation, causing the highly conserved DFG-motif (Asp(1150), Phe(1151), and Gly(1152)) to switch from the inactive "D-out/F-in" to the nucleotide-binding-competent "D-in/F-out" conformation. The minimum free-energy path computed from the string method preserves these helical intermediates along the inactive-to-active path, and the thermodynamic free-energy differences are consistent with previous work on various other kinases. The mechanisms revealed by TAMD also suggest that the regulatory spine can be dynamically assembled/disassembled either by DFG-flip or by movement of the αC-helix. Together, these findings both broaden our understanding of kinase activation and point to intermediates as specific therapeutic targets.

Comparison of Three Chain-of-States Methods: Nudged Elastic Band and Replica Path with Restraints or Constraints

Chain-of-state methods are becoming important tools in studying the chemical reaction mechanisms, especially for biomacromolecules. In this article, three chain-of-state methods, nudged elastic band (NEB) method and the replica path method with restraints or constraints, were tested and compared using three model systems with various sizes and at different levels of theory: alanine dipeptide isomerization, β-alanine intramolecular condensation, and the matrix metalloproteinase 2 inhibition mechanism. The levels of theory used to describe the three model systems include molecular mechanics (MM), quantum mechanics (QM), and combined quantum mechanics and molecular mechanics (QM/MM). All three methods could correctly determine a reaction path with reasonable estimation of reaction barriers in most cases. The RMSD measurement with additional weighting schemes provides practically infinite choices of reaction coordinates to describe the reaction progress. These findings demonstrate that the chain-of-state methods are powerful tools when being used carefully to generate a plausible reaction mechanism with full pathway for complex systems at an affordable computational cost.

Unrestrained computation of free energy along a path

We apply the adaptive biasing potential (ABP) method to optimize the principal curve defining a conformational transition between two known end states and to subsequently compute the one-dimensional potential of mean force as a function of arc length along the principal curve. This approach allows the use of the ABP method in a collective variable space of arbitrary dimension and offers several advantages over line-search methods. First, configurations are neither generated along an initial path for the transition nor equilibrated during evolution of the path. Second, and most importantly, the powerful sampling provided by the ABP serves to accelerate the dynamics during the optimization and computation of the free energy. Finally, the free energy is formulated as a potential of mean force that captures changes in the reaction channel along the principal curve, in contrast to the free energy profile evaluated from the local free-energy gradient in restrained path optimization methods. We first demonstrate the ABP formulation of path optimization using a two-dimensional potential surface and then with a more complex system of Src protein tyrosine kinase. The method is shown to be efficient and robust in the case of rugged, free-energy landscapes.

Energetics and mechanism of the normal-to-amyloidogenic isomerization of β2-microglobulin: on-the-fly string method calculations

We use on-the-fly finite temperature string method in collective variables to study the transition from a normal to an amyloidogenic conformation of β2-microglobulin. We show that the protonation state of two histidine residues is of key importance and that under acidic (protonating) conditions, the transition to the amyloidgenic form is facilitated by both displacement of N-terminal residues to disrupt a hydrophobic pocket and by side-chain/side-chain electrostatic attraction, both of which facilitate a cis-trans prolyl isomerization. The free energy barriers for the normal-to-amyloidogenic isomerization are found to be 14.9 and 7.1 kcal/mol for the neutral and protonated cases, respectively, consistent with enhanced amyloidgenesis at low pH observed both in vitro and in hemodialysis-associated amyloidosis, and somewhat lower than experimentally determined barriers for bare prolyl cis-trans isomerization. We suggest specific mutagenesis experiments which could be used to further validate the mechanism observed.