Reconstruction of the Most Probable Folding Transition Path from All-Atom Replica Exchange Simulations, using the Dynamic String Method

Camilo A. Jimenez-Cruz, Angel E. Garcia
2013 Biophysical Journal  
formed, cyt c becomes a peroxidase which catalyzes peroxidation of polyunsaturated CL species into products that are required for the execution of the apoptotic program. The essential details of the cyt c/CL interactions -leading to the gain of cyt c's peroxidase competence towards CL -remain to be elucidated. Here, we used Coarse Graining Molecular Dynamic (CG-MD) simulation as a computational route to explore structural details of cyt c interactions with CL-containing phospholipid bilayer.
more » ... ng the MARTINI force field of a 4-to-1 atoms-to-bead mapping for both lipids and protein, we simulated a lipid bilayer containing 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) and 20% tetraoleoyl-CL (TOCL). Furthermore, we investigated the effect of CL peroxidation on membrane structure and interactions with cyt c. The simulation visualizes cyt c movement toward the phospholipid bilayer, to which it binds and then migrates along its surface. Inside the bilayer, clustering of CL with respect to the position of cyt c on the bilayer surface was observed. Peroxidation of CL altered the phospholipid membrane and caused the appearance of disorganized domains. Furthermore, increasing levels of peroxidation resulted in decreasing capacity of cyt c to interact with the membrane, whereby cyt c completely lost the ability to bind to the bilayer containing fully peroxidized CL (with all four of CL's acyls oxidiatively modified). Large-scale conformational transitions in biomolecules occur on complex, high-dimensional free energy surfaces. Simulating these transitions is often challenging due to the mismatch between the timescales required to traverse the energetic landscape and those accessible by conventional simulation techniques. The "weighted ensemble" (WE) path sampling method is a rigorous technique for simulating this class of rare transitions. We extend the WE method by combining it with a string method to adaptively refine the set of order parameters used to enhance sampling along the transition pathway. This permits sampling transitions in the space of many order parameters for a wide range of equilibrium and non-equilibrium processes. From these simulations, accurate estimates of steady state conformational distributions and reaction rates can be obtained, even for systems with complex transition pathways that may involve metastable intermediates. We demonstrate the application of this method using two simple models of driven Brownian motion, protein conformational change with a coarse-grained model, and fully-atomistic models. Due to the challenges involved with modeling complex molecular systems,it is essential that computationally intelligent schemes be produced which put the computational effort where and when it is needed to capture important phenomena, and maintain needed accuracy at minimum costs. In this work, we develop and investigate algorithms for the adaptive modeling and simulation of the dynamic behavior of highly complex multiscale processes. This is accomplished through the appropriate use of an adaptive hybridization of existing, newly developed, and proposed advanced multibody dynamics algorithms and modeling strategies for forward dynamic simulation. The adaptive multiscale simulation technique presented here benefits from the highly parallelizable structure of the divide and conquer (DCA) framework for modeling multibody systems. These algorithms permits a large complex molecule (or systems of molecules) to be seamlessly treated using a hierarchy of reduced order models ranging from atomistic to the continuum scale. The reduced order and low fidelity models, when correctly developed, can provide significant computational savings. The reason is these coarser scale models effectively constrain out high frequency modes which dominate the integration of the equations of motion, but are of little or no relevance to the important overall conformational behavior of the system. When such fine-scale (temporal and/or spatial) scales are not needed the subassemblies of the molecule can be identified and modeled by coarser scale representations which may include rigidbody models, articulated body models, flexible body models, and continuum models. The simulation technique using the DCA framework permits switching between different resolution models adaptively during the simulation. For example, from fully atomistic to a multi-flexible or multi-rigid body representation can be achieved depending on the researcher specified internal metric indicators or error estimates. We have developed LOOS (Lightweight Object Oriented Structure-analysis) as a tool for making new tools for analyzing molecular simulations. LOOS is an object-oriented library designed to facilitate the rapid development of new methods for structural analysis. LOOS includes over 120 pre-built tools for common structural analysis tasks, such as assessing the convergence of simulations, hydrogen bonding patterns, and the construction and visualization of 3D histograms as density. LOOS supports reading the native file formats of most common simulation packages, including AMBER, CHARMM, Gromacs, NAMD, and Tinker and can write NAMD formats (PDB and DCD). A dynamic atom selection language, based on the C expression syntax, is included as part of the library and is easily accessible to both the programmer and the end user. LOOS is written in Cþþ and makes extensive use of Boost and the Standard Template Library. Through modern Cþþ design, LOOS is both simple to use (requiring knowledge of only 4 core classes and a few utility functions) and easily extensible. A Python interface to the core components of LOOS is also available, further facilitating rapid development of analysis tools and broadening the LOOS community by making it accessible to those who would otherwise be deterred by using Cþþ. LOOS also includes a set of libraries and tools for performing elastic network model calculations that are easily extended to accommodate new methods. In this work we demonstrate how the most probable transition path between metastable states can be recovered from replica exchange molecular dynamic simulation data by using the dynamic string method. The local drift vector in collective variables is determined via the short continuous trajectories between replica exchanges at a given temperature. A string is updated based on this drift vector to produce reaction pathways between the folded and unfolded state. The method is applied to a beta hairpin-forming peptide to obtain information on the folding mechanism and transition state using different sets of collective variables. Two main folding pathways differing in the order of events are found and discussed, and the relative free energy differences for each path estimated. The string that shows turn formation before native arrangement of the hydrophobic core has lower free energy difference from the unfolded state and thus is expected to be more probable. Finally, in both cases, the conformations on the peak of the free energy barrier have the tryptophans on the same side of the beta structure suggesting this arrangement to be the rate-limiting step. The ensemble of 3-d configurations exhibited by a biomolecule, that is, its intrinsic motion, can be altered by several environmental factors, and also by the binding of other molecules. Quantification of such induced changes in intrinsic motion is important because it provides a basis for relating thermodynamic changes to changes in molecular motion. This task is, however, challenging because it requires comparing two high-dimensional datasets. Traditionally, when analyzing molecular simulations, this problem is circumvented by first reducing the dimensions of the two ensembles separately, and then comparing summary statistics from the two ensembles against each other. However, since dimensionality reduction is carried out prior to ensemble comparison, such strategies are susceptible to artifactual biases from information loss. Here we introduce a method based on support vector machines that yields a normalized quantitative estimate for the difference between two ensembles after comparing them directly against one another. While this method can be applied to any molecular system, including non-biological molecules and crystals, here we show how it can be applied to identify the specific regions of a paramyxovirus G-protein that are affected by the binding of its preferred human receptor, Ephrin-B2. This protein-protein interaction essentially regulates viral fusion. Specifically, for every residue in the G-protein, we obtain separately a quantitative 504a
doi:10.1016/j.bpj.2012.11.2784 fatcat:zicgixufkzduril6zqd4sqkhpe