1
|
Sun Y. Efficient acceleration of the convergence of the minimum free energy path via a path-planning generated initial guess. J Comput Chem 2024. [PMID: 39291721 DOI: 10.1002/jcc.27504] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 05/17/2024] [Revised: 07/25/2024] [Accepted: 08/29/2024] [Indexed: 09/19/2024]
Abstract
We demonstrate that combining a shifted clustering algorithm with a fast-marching-based algorithm can generate accurate approximations of the minimum energy path (MEP) given a free energy landscape (FEL). Using this approximation as the initial guess for the MEP, followed by further refinement with the string method (referred to as the fast marching tree (FMT)-string combined approach), significantly reduces the number of iterations required for MEP convergence. This approach saves substantial time compared to using linear interpolation (LI) for the initial guess. Our method offers a viable solution for obtaining an effective initial guess of the MEP when an approximate or converged FEL is available. This work highlights the potential of applying FMT-based approaches to extract the MEP in chemical reactions.
Collapse
Affiliation(s)
- Yi Sun
- Department of Chemistry, Chicago Center for Theoretical Chemistry, James Franck Institute, and Institute for Biophysical Dynamics, The University of Chicago, Chicago, Illinois, USA
| |
Collapse
|
2
|
Jiang W. Studying the Collective Functional Response of a Receptor in Alchemical Ligand Binding Free Energy Simulations with Accelerated Solvation Layer Dynamics. J Chem Theory Comput 2024; 20:3085-3095. [PMID: 38568961 DOI: 10.1021/acs.jctc.4c00191] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [MESH Headings] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 04/05/2024]
Abstract
Ligand binding free energy simulations (LB-FES) that involve sampling of protein functional conformations have been longstanding challenges in research on molecular recognition. Particularly, modeling of the conformational transition pathway and design of the heuristic biasing mechanism are severe bottlenecks for the existing enhanced configurational sampling (ECS) methods. Inspired by the key role of hydration in regulating conformational dynamics of macromolecules, this report proposes a novel ECS approach that facilitates binding-associated structural dynamics by accelerated hydration transitions in combination with the λ-exchange of free energy perturbation (FEP). Two challenging protein-ligand binding processes involving large configurational transitions of the receptor are studied, with hydration transitions at binding sites accelerated by Hamiltonian-simulated annealing of the hydration layer. Without the need for pathway analysis or ad hoc barrier flattening potential, LB-FES were performed with FEP/λ-exchange molecular dynamics simulation at a minor overhead for annealing of the hydration layer. The LB-FES studies showed that the accelerated rehydration significantly enhances the collective conformational transitions of the receptor, and convergence of binding affinity calculations is obtained at a sweet-spot simulation time scale. Alchemical LB-FES with the proposed ECS strategy is free from the effort of trial and error for the setup and realizes efficient on-the-fly sampling for the collective functional response of the receptor and bound water and therefore presents a practical approach to high-throughput screening in drug discovery.
Collapse
Affiliation(s)
- Wei Jiang
- Computational Science Division, Argonne National Laboratory, 9700 South Cass Avenue, Building 240, Argonne, Illinois 60439, United States
| |
Collapse
|
3
|
Blazhynska M, Gumbart JC, Chen H, Tajkhorshid E, Roux B, Chipot C. A Rigorous Framework for Calculating Protein-Protein Binding Affinities in Membranes. J Chem Theory Comput 2023; 19:9077-9092. [PMID: 38091976 PMCID: PMC11145395 DOI: 10.1021/acs.jctc.3c00941] [Citation(s) in RCA: 2] [Impact Index Per Article: 2.0] [Reference Citation Analysis] [Abstract] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/27/2023]
Abstract
Calculating the binding free energy of integral transmembrane (TM) proteins is crucial for understanding the mechanisms by which they recognize one another and reversibly associate. The glycophorin A (GpA) homodimer, composed of two α-helical segments, has long served as a model system for studying TM protein reversible association. The present work establishes a methodological framework for calculating the binding affinity of the GpA homodimer in the heterogeneous environment of a membrane. Our investigation carefully considered a variety of protocols, including the appropriate choice of the force field, rigorous standardization reflecting the experimental conditions, sampling algorithm, anisotropic environment, and collective variables, to accurately describe GpA dimerization via molecular dynamics-based approaches. Specifically, two strategies were explored: (i) an unrestrained potential mean force (PMF) calculation, which merely enhances sampling along the separation of the two binding partners without any restraint, and (ii) a so-called "geometrical route", whereby the α-helices are progressively separated with imposed restraints on their orientational, positional, and conformational degrees of freedom to accelerate convergence. Our simulations reveal that the simplified, unrestrained PMF approach is inadequate for the description of GpA dimerization. Instead, the geometrical route, tailored specifically to GpA in a membrane environment, yields excellent agreement with experimental data within a reasonable computational time. A dimerization free energy of -10.7 kcal/mol is obtained, in fairly good agreement with available experimental data. The geometrical route further helps elucidate how environmental forces drive association before helical interactions stabilize it. Our simulations also brought to light a distinct, long-lived spatial arrangement that potentially serves as an intermediate state during dimer formation. The methodological advances in the generalized geometrical route provide a powerful tool for accurate and efficient binding-affinity calculations of intricate TM protein complexes in inhomogeneous environments.
Collapse
Affiliation(s)
- Marharyta Blazhynska
- Laboratoire International Associé Centre National de la Recherche Scientifique et University of Illinois at Urbana-Champaign, Unité Mixte de Recherche n°7019, Université de Lorraine, B.P. 70239, Vandœuvre-lès-Nancy cedex 54506, France
| | - James C Gumbart
- School of Physics, Georgia Institute of Technology, 837 State Street, Atlanta, Georgia 30332, United States
| | - Haochuan Chen
- Laboratoire International Associé Centre National de la Recherche Scientifique et University of Illinois at Urbana-Champaign, Unité Mixte de Recherche n°7019, Université de Lorraine, B.P. 70239, Vandœuvre-lès-Nancy cedex 54506, France
| | - Emad Tajkhorshid
- Theoretical and Computational Biophysics Group, NIH Center for Macromolecular Modeling and Visualization, Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, 405 N. Mathews Avenue, Urbana, Illinois 61801, United States
- Department of Biochemistry, University of Illinois at Urbana-Champaign, 600 S. Mathews Avenue, Urbana, Illinois 61801, United States
| | - Benoît Roux
- Department of Biochemistry and Molecular Biology, The University of Chicago, 929 E. 57th Street W225, Chicago, Illinois 60637, United States
- Department of Chemistry, The University of Chicago, 5735 S. Ellis Avenue, Chicago, Illinois 60637, United States
| | - Christophe Chipot
- Laboratoire International Associé Centre National de la Recherche Scientifique et University of Illinois at Urbana-Champaign, Unité Mixte de Recherche n°7019, Université de Lorraine, B.P. 70239, Vandœuvre-lès-Nancy cedex 54506, France
- Theoretical and Computational Biophysics Group, NIH Center for Macromolecular Modeling and Visualization, Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, 405 N. Mathews Avenue, Urbana, Illinois 61801, United States
- Department of Biochemistry and Molecular Biology, The University of Chicago, 929 E. 57th Street W225, Chicago, Illinois 60637, United States
- Department of Chemistry, The University of Hawai'i at Ma̅noa, 2545 McCarthy Mall, Honolulu, Hawaii 96822, United States
| |
Collapse
|
4
|
Ojha AA, Votapka LW, Amaro RE. QMrebind: incorporating quantum mechanical force field reparameterization at the ligand binding site for improved drug-target kinetics through milestoning simulations. Chem Sci 2023; 14:13159-13175. [PMID: 38023523 PMCID: PMC10664576 DOI: 10.1039/d3sc04195f] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Grants] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 08/11/2023] [Accepted: 10/22/2023] [Indexed: 12/01/2023] Open
Abstract
Understanding the interaction of ligands with biomolecules is an integral component of drug discovery and development. Challenges for computing thermodynamic and kinetic quantities for pharmaceutically relevant receptor-ligand complexes include the size and flexibility of the ligands, large-scale conformational rearrangements of the receptor, accurate force field parameters, simulation efficiency, and sufficient sampling associated with rare events. Our recently developed multiscale milestoning simulation approach, SEEKR2 (Simulation Enabled Estimation of Kinetic Rates v.2), has demonstrated success in predicting unbinding (koff) kinetics by employing molecular dynamics (MD) simulations in regions closer to the binding site. The MD region is further subdivided into smaller Voronoi tessellations to improve the simulation efficiency and parallelization. To date, all MD simulations are run using general molecular mechanics (MM) force fields. The accuracy of calculations can be further improved by incorporating quantum mechanical (QM) methods into generating system-specific force fields through reparameterizing ligand partial charges in the bound state. The force field reparameterization process modifies the potential energy landscape of the bimolecular complex, enabling a more accurate representation of the intermolecular interactions and polarization effects at the bound state. We present QMrebind (Quantum Mechanical force field reparameterization at the receptor-ligand binding site), an ORCA-based software that facilitates reparameterizing the potential energy function within the phase space representing the bound state in a receptor-ligand complex. With SEEKR2 koff estimates and experimentally determined kinetic rates, we compare and interpret the receptor-ligand unbinding kinetics obtained using the newly reparameterized force fields for model host-guest systems and HSP90-inhibitor complexes. This method provides an opportunity to achieve higher accuracy in predicting receptor-ligand koff rate constants.
Collapse
Affiliation(s)
- Anupam Anand Ojha
- Department of Chemistry and Biochemistry, University of California San Diego La Jolla California 92093 USA
| | - Lane William Votapka
- Department of Chemistry and Biochemistry, University of California San Diego La Jolla California 92093 USA
| | - Rommie Elizabeth Amaro
- Department of Molecular Biology, University of California San Diego La Jolla California 92093 USA
| |
Collapse
|
5
|
Rizzi V, Aureli S, Ansari N, Gervasio FL. OneOPES, a Combined Enhanced Sampling Method to Rule Them All. J Chem Theory Comput 2023; 19:5731-5742. [PMID: 37603295 PMCID: PMC10500989 DOI: 10.1021/acs.jctc.3c00254] [Citation(s) in RCA: 1] [Impact Index Per Article: 1.0] [Reference Citation Analysis] [Abstract] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 03/06/2023] [Indexed: 08/22/2023]
Abstract
Enhanced sampling techniques have revolutionized molecular dynamics (MD) simulations, enabling the study of rare events and the calculation of free energy differences in complex systems. One of the main families of enhanced sampling techniques uses physical degrees of freedom called collective variables (CVs) to accelerate a system's dynamics and recover the original system's statistics. However, encoding all the relevant degrees of freedom in a limited number of CVs is challenging, particularly in large biophysical systems. Another category of techniques, such as parallel tempering, simulates multiple replicas of the system in parallel, without requiring CVs. However, these methods may explore less relevant high-energy portions of the phase space and become computationally expensive for large systems. To overcome the limitations of both approaches, we propose a replica exchange method called OneOPES that combines the power of multireplica simulations and CV-based enhanced sampling. This method efficiently accelerates the phase space sampling without the need for ideal CVs, extensive parameters fine tuning nor the use of a large number of replicas, as demonstrated by its successful applications to protein-ligand binding and protein folding benchmark systems. Our approach shows promise as a new direction in the development of enhanced sampling techniques for molecular dynamics simulations, providing an efficient and robust framework for the study of complex and unexplored problems.
Collapse
Affiliation(s)
- Valerio Rizzi
- School
of Pharmaceutical Sciences, University of
Geneva, Rue Michel Servet 1, 1206 Genève, Switzerland
- Institute
of Pharmaceutical Sciences of Western Switzerland (ISPSO), University of Geneva, 1206 Genève, Switzerland
- Swiss
Institute of Bioinformatics, University
of Geneva, 1206 Genève, Switzerland
| | - Simone Aureli
- School
of Pharmaceutical Sciences, University of
Geneva, Rue Michel Servet 1, 1206 Genève, Switzerland
- Institute
of Pharmaceutical Sciences of Western Switzerland (ISPSO), University of Geneva, 1206 Genève, Switzerland
- Swiss
Institute of Bioinformatics, University
of Geneva, 1206 Genève, Switzerland
| | - Narjes Ansari
- Atomistic
Simulations, Italian Institute of Technology, Via Enrico Melen 83, 16152 Genova, Italy
| | - Francesco Luigi Gervasio
- School
of Pharmaceutical Sciences, University of
Geneva, Rue Michel Servet 1, 1206 Genève, Switzerland
- Institute
of Pharmaceutical Sciences of Western Switzerland (ISPSO), University of Geneva, 1206 Genève, Switzerland
- Swiss
Institute of Bioinformatics, University
of Geneva, 1206 Genève, Switzerland
- Department
of Chemistry, University College London, WC1E 6BT London, U.K.
| |
Collapse
|
6
|
Saurabh S, Nadendla K, Purohit SS, Sivakumar PM, Cetinel S. Fuzzy Drug Targets: Disordered Proteins in the Drug-Discovery Realm. ACS OMEGA 2023; 8:9729-9747. [PMID: 36969402 PMCID: PMC10034788 DOI: 10.1021/acsomega.2c07708] [Citation(s) in RCA: 8] [Impact Index Per Article: 8.0] [Reference Citation Analysis] [Abstract] [Grants] [Track Full Text] [Figures] [Subscribe] [Scholar Register] [Received: 12/02/2022] [Accepted: 02/17/2023] [Indexed: 06/18/2023]
Abstract
Intrinsically disordered proteins (IDPs) and regions (IDRs) form a large part of the eukaryotic proteome. Contrary to the structure-function paradigm, the disordered proteins perform a myriad of functions in vivo. Consequently, they are involved in various disease pathways and are plausible drug targets. Unlike folded proteins, that have a defined structure and well carved out drug-binding pockets that can guide lead molecule selection, the disordered proteins require alternative drug-development methodologies that are based on an acceptable picture of their conformational ensemble. In this review, we discuss various experimental and computational techniques that contribute toward understanding IDP "structure" and describe representative pursuances toward IDP-targeting drug development. We also discuss ideas on developing rational drug design protocols targeting IDPs.
Collapse
Affiliation(s)
- Suman Saurabh
- Molecular
Sciences Research Hub, Department of Chemistry, Imperial College London, London W12 0BZ, U.K.
| | - Karthik Nadendla
- Center
for Misfolding Diseases, Yusuf Hamied Department of Chemistry, Lensfield
Road, University of Cambridge, Cambridge CB2 1EW, U.K.
| | - Shubh Sanket Purohit
- Department
of Clinical Haematology, Sahyadri Superspeciality
Hospital, Pune, Maharashtra 411038, India
| | - Ponnurengam Malliappan Sivakumar
- Institute
of Research and Development, Duy Tan University, Da Nang 550000, Vietnam
- School
of Medicine and Pharmacy, Duy Tan University, Da Nang 550000, Vietnam
- Nanotechnology
Research and Application Center (SUNUM), Sabanci University, Istanbul 34956, Turkey
| | - Sibel Cetinel
- Nanotechnology
Research and Application Center (SUNUM), Sabanci University, Istanbul 34956, Turkey
- Faculty of
Engineering and Natural Sciences, Molecular Biology, Genetics and
Bioengineering Program, Sabanci University, Istanbul 34956, Turkey
| |
Collapse
|
7
|
Koskin V, Kells A, Clayton J, Hartmann AK, Annibale A, Rosta E. Variational kinetic clustering of complex networks. J Chem Phys 2023; 158:104112. [PMID: 36922127 DOI: 10.1063/5.0105099] [Citation(s) in RCA: 1] [Impact Index Per Article: 1.0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/14/2022] Open
Abstract
Efficiently identifying the most important communities and key transition nodes in weighted and unweighted networks is a prevalent problem in a wide range of disciplines. Here, we focus on the optimal clustering using variational kinetic parameters, linked to Markov processes defined on the underlying networks, namely, the slowest relaxation time and the Kemeny constant. We derive novel relations in terms of mean first passage times for optimizing clustering via the Kemeny constant and show that the optimal clustering boundaries have equal round-trip times to the clusters they separate. We also propose an efficient method that first projects the network nodes onto a 1D reaction coordinate and subsequently performs a variational boundary search using a parallel tempering algorithm, where the variational kinetic parameters act as an energy function to be extremized. We find that maximization of the Kemeny constant is effective in detecting communities, while the slowest relaxation time allows for detection of transition nodes. We demonstrate the validity of our method on several test systems, including synthetic networks generated from the stochastic block model and real world networks (Santa Fe Institute collaboration network, a network of co-purchased political books, and a street network of multiple cities in Luxembourg). Our approach is compared with existing clustering algorithms based on modularity and the robust Perron cluster analysis, and the identified transition nodes are compared with different notions of node centrality.
Collapse
Affiliation(s)
- Vladimir Koskin
- Department of Chemistry, King's College London, SE1 1DB London, United Kingdom
| | - Adam Kells
- Department of Chemistry, King's College London, SE1 1DB London, United Kingdom
| | - Joe Clayton
- Department of Physics and Astronomy, University College London, WC1E 6BT London, United Kingdom
| | | | - Alessia Annibale
- Department of Mathematics, King's College London, SE11 6NJ London, United Kingdom
| | - Edina Rosta
- Department of Physics and Astronomy, University College London, WC1E 6BT London, United Kingdom
| |
Collapse
|
8
|
Buslaev P, Aho N, Jansen A, Bauer P, Hess B, Groenhof G. Best Practices in Constant pH MD Simulations: Accuracy and Sampling. J Chem Theory Comput 2022; 18:6134-6147. [PMID: 36107791 PMCID: PMC9558372 DOI: 10.1021/acs.jctc.2c00517] [Citation(s) in RCA: 9] [Impact Index Per Article: 4.5] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/11/2023]
Abstract
![]()
Various approaches
have been proposed to include the
effect of
pH in molecular dynamics (MD) simulations. Among these, the λ-dynamics approach proposed
by Brooks and
co-workers [Kong, X.; Brooks III, C. L. J. Chem. Phys.1996, 105, 2414−2423] can be performed
with little computational overhead and hfor each typeence be used
to routinely perform MD simulations at microsecond time scales, as
shown in the accompanying paper [Aho, N. et al. J. Chem. Theory
Comput.2022, DOI: 10.1021/acs.jctc.2c00516]. At
such time scales, however, the accuracy of the molecular mechanics
force field and the parametrization becomes critical. Here, we address
these issues and provide the community with guidelines on how to set
up and perform long time scale constant pH MD simulations. We found
that barriers associated with the torsions of side chains in the CHARMM36m
force field are too high for reaching convergence in constant pH MD
simulations on microsecond time scales. To avoid the high computational
cost of extending the sampling, we propose small modifications to
the force field to selectively reduce the torsional barriers. We demonstrate
that with such modifications we obtain converged distributions of
both protonation and torsional degrees of freedom and hence consistent
pKa estimates, while the sampling of the
overall configurational space accessible to proteins is unaffected
as compared to normal MD simulations. We also show that the results
of constant pH MD depend on the accuracy of the correction potentials.
While these potentials are typically obtained by fitting a low-order
polynomial to calculated free energy profiles, we find that higher
order fits are essential to provide accurate and consistent results.
By resolving problems in accuracy and sampling, the work described
in this and the accompanying paper paves the way to the widespread
application of constant pH MD beyond pKa prediction.
Collapse
Affiliation(s)
- Pavel Buslaev
- Nanoscience Center and Department of Chemistry, University of Jyväskylä, 40014 Jyväskylä, Finland
| | - Noora Aho
- Nanoscience Center and Department of Chemistry, University of Jyväskylä, 40014 Jyväskylä, Finland
| | - Anton Jansen
- Department of Applied Physics, Science for Life Laboratory, KTH Royal Institute of Technology, 100 44 Stockholm, Sweden
| | - Paul Bauer
- Department of Applied Physics, Science for Life Laboratory, KTH Royal Institute of Technology, 100 44 Stockholm, Sweden
| | - Berk Hess
- Department of Applied Physics and Swedish e-Science Research Center, Science for Life Laboratory, KTH Royal Institute of Technology, 100 44 Stockholm, Sweden
| | - Gerrit Groenhof
- Nanoscience Center and Department of Chemistry, University of Jyväskylä, 40014 Jyväskylä, Finland
| |
Collapse
|
9
|
Modeling Adsorption, Conformation, and Orientation of the Fis1 Tail Anchor at the Mitochondrial Outer Membrane. MEMBRANES 2022; 12:membranes12080752. [PMID: 36005667 PMCID: PMC9413518 DOI: 10.3390/membranes12080752] [Citation(s) in RCA: 1] [Impact Index Per Article: 0.5] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Subscribe] [Scholar Register] [Received: 06/29/2022] [Revised: 07/28/2022] [Accepted: 07/28/2022] [Indexed: 12/04/2022]
Abstract
Proteins can be targeted to organellar membranes by using a tail anchor (TA), a stretch of hydrophobic amino acids found at the polypeptide carboxyl-terminus. The Fis1 protein (Fis1p), which promotes mitochondrial and peroxisomal division in the yeast Saccharomyces cerevisiae, is targeted to those organelles by its TA. Substantial evidence suggests that Fis1p insertion into the mitochondrial outer membrane can occur without the need for a translocation machinery. However, recent findings raise the possibility that Fis1p insertion into mitochondria might be promoted by a proteinaceous complex. Here, we have performed atomistic and coarse-grained molecular dynamics simulations to analyze the adsorption, conformation, and orientation of the Fis1(TA). Our results support stable insertion at the mitochondrial outer membrane in a monotopic, rather than a bitopic (transmembrane), configuration. Once inserted in the monotopic orientation, unassisted transition to the bitopic orientation is expected to be blocked by the highly charged nature of the TA carboxyl-terminus and by the Fis1p cytosolic domain. Our results are consistent with a model in which Fis1p does not require a translocation machinery for insertion at mitochondria.
Collapse
|
10
|
Raubenolt BA, Islam NN, Summa CM, Rick SW. Molecular dynamics simulations of the flexibility and inhibition of SARS-CoV-2 NSP 13 helicase. J Mol Graph Model 2022; 112:108122. [PMID: 35021142 PMCID: PMC8730789 DOI: 10.1016/j.jmgm.2022.108122] [Citation(s) in RCA: 3] [Impact Index Per Article: 1.5] [Reference Citation Analysis] [Abstract] [MESH Headings] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 09/30/2021] [Revised: 12/31/2021] [Accepted: 01/03/2022] [Indexed: 11/25/2022]
Abstract
The helicase protein of the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) is both a good potential drug target and very flexible. The flexibility, and therefore its function, could be reduced through knowledge of these motions and identification of allosteric pockets. Using molecular dynamics simulations with enhanced sampling, we determined key modes of motion and sites on the protein that are at the interface between flexible domains of the proteins. We developed an approach to map the principal components of motion onto the surface of a potential binding pocket to help in the identification of allosteric sites.
Collapse
Affiliation(s)
- Bryan A Raubenolt
- Department of Chemistry, University of New Orleans, New Orleans, LA, 70148, USA
| | - Naeyma N Islam
- Department of Chemistry, University of New Orleans, New Orleans, LA, 70148, USA
| | - Christoper M Summa
- Department of Computer Science, University of New Orleans, New Orleans, LA, 70148, USA
| | - Steven W Rick
- Department of Chemistry, University of New Orleans, New Orleans, LA, 70148, USA.
| |
Collapse
|
11
|
Xie H, Rojas A, Maisuradze GG, Khelashvili G. Mechanistic Kinetic Model Reveals How Amyloidogenic Hydrophobic Patches Facilitate the Amyloid-β Fibril Elongation. ACS Chem Neurosci 2022; 13:987-1001. [PMID: 35258946 PMCID: PMC8986627 DOI: 10.1021/acschemneuro.1c00801] [Citation(s) in RCA: 3] [Impact Index Per Article: 1.5] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/29/2022] Open
Abstract
Abnormal aggregation of amyloid β (Aβ) peptides into fibrils plays a critical role in the development of Alzheimer's disease. A two-stage "dock-lock" model has been proposed for the Aβ fibril elongation process. However, the mechanisms of the Aβ monomer-fibril binding process have not been elucidated with the necessary molecular-level precision, so it remains unclear how the lock phase dynamics leads to the overall in-register binding of the Aβ monomer onto the fibril. To gain mechanistic insights into this critical step during the fibril elongation process, we used molecular dynamics (MD) simulations with a physics-based coarse-grained UNited-RESidue (UNRES) force field and sampled extensively the dynamics of the lock phase process, in which a fibril-bound Aβ(9-40) peptide rearranged to establish the native docking conformation. Analysis of the MD trajectories with Markov state models was used to quantify the kinetics of the rearrangement process and the most probable pathways leading to the overall native docking conformation of the incoming peptide. These revealed a key intermediate state in which an intra-monomer hairpin is formed between the central core amyloidogenic patch 18VFFA21 and the C-terminal hydrophobic patch 34LMVG37. This hairpin structure is highly favored as a transition state during the lock phase of the fibril elongation. We propose a molecular mechanism for facilitation of the Aβ fibril elongation by amyloidogenic hydrophobic patches.
Collapse
Affiliation(s)
- Hengyi Xie
- Department of Physiology and Biophysics, Weill Cornell Medicine, New York, New York 10065, United States
| | - Ana Rojas
- Schrödinger, Inc., 1540 Broadway, 24th Floor, New York, New York 10036, United States
| | - Gia G. Maisuradze
- Baker Laboratory of Chemistry and Chemical Biology, Cornell University, Ithaca, New York 14853, United States
| | - George Khelashvili
- Department of Physiology and Biophysics, Weill Cornell Medicine, New York, New York 10065, United States
| |
Collapse
|
12
|
Mlýnský V, Janeček M, Kührová P, Fröhlking T, Otyepka M, Bussi G, Banáš P, Šponer J. Toward Convergence in Folding Simulations of RNA Tetraloops: Comparison of Enhanced Sampling Techniques and Effects of Force Field Modifications. J Chem Theory Comput 2022; 18:2642-2656. [PMID: 35363478 DOI: 10.1021/acs.jctc.1c01222] [Citation(s) in RCA: 17] [Impact Index Per Article: 8.5] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/29/2022]
Abstract
Atomistic molecular dynamics simulations represent an established technique for investigation of RNA structural dynamics. Despite continuous development, contemporary RNA simulations still suffer from suboptimal accuracy of empirical potentials (force fields, ffs) and sampling limitations. Development of efficient enhanced sampling techniques is important for two reasons. First, they allow us to overcome the sampling limitations, and second, they can be used to quantify ff imbalances provided they reach a sufficient convergence. Here, we study two RNA tetraloops (TLs), namely the GAGA and UUCG motifs. We perform extensive folding simulations and calculate folding free energies (ΔGfold°) with the aim to compare different enhanced sampling techniques and to test several modifications of the nonbonded terms extending the AMBER OL3 RNA ff. We demonstrate that replica-exchange solute tempering (REST2) simulations with 12-16 replicas do not show any sign of convergence even when extended to a timescale of 120 μs per replica. However, the combination of REST2 with well-tempered metadynamics (ST-MetaD) achieves good convergence on a timescale of 5-10 μs per replica, improving the sampling efficiency by at least 2 orders of magnitude. Effects of ff modifications on ΔGfold° energies were initially explored by the reweighting approach and then validated by new simulations. We tested several manually prepared variants of the gHBfix potential which improve stability of the native state of both TLs by ∼2 kcal/mol. This is sufficient to conveniently stabilize the folded GAGA TL while the UUCG TL still remains under-stabilized. Appropriate adjustment of van der Waals parameters for C-H···O5' base-phosphate interaction may further stabilize the native states of both TLs by ∼0.6 kcal/mol.
Collapse
Affiliation(s)
- Vojtěch Mlýnský
- Institute of Biophysics of the Czech Academy of Sciences, Královopolská 135, 612 65 Brno, Czech Republic
| | - Michal Janeček
- Department of Physical Chemistry, Faculty of Science, Palacký University, tř. 17 listopadu 12, 771 46 Olomouc, Czech Republic
| | - Petra Kührová
- Institute of Biophysics of the Czech Academy of Sciences, Královopolská 135, 612 65 Brno, Czech Republic.,Regional Centre of Advanced Technologies and Materials, Czech Advanced Technology and Research Institute (CATRIN), Palacký University Olomouc, Šlechtitelů 27, 779 00 Olomouc, Czech Republic
| | - Thorben Fröhlking
- Scuola Internazionale Superiore di Studi Avanzati, SISSA, via Bonomea 265, 34136 Trieste, Italy
| | - Michal Otyepka
- Regional Centre of Advanced Technologies and Materials, Czech Advanced Technology and Research Institute (CATRIN), Palacký University Olomouc, Šlechtitelů 27, 779 00 Olomouc, Czech Republic.,IT4Innovations, VSB─Technical University of Ostrava, 17. listopadu 2172/15, 708 00 Ostrava-Poruba, Czech Republic
| | - Giovanni Bussi
- Scuola Internazionale Superiore di Studi Avanzati, SISSA, via Bonomea 265, 34136 Trieste, Italy
| | - Pavel Banáš
- Institute of Biophysics of the Czech Academy of Sciences, Královopolská 135, 612 65 Brno, Czech Republic.,Regional Centre of Advanced Technologies and Materials, Czech Advanced Technology and Research Institute (CATRIN), Palacký University Olomouc, Šlechtitelů 27, 779 00 Olomouc, Czech Republic
| | - Jiří Šponer
- Institute of Biophysics of the Czech Academy of Sciences, Královopolská 135, 612 65 Brno, Czech Republic.,Regional Centre of Advanced Technologies and Materials, Czech Advanced Technology and Research Institute (CATRIN), Palacký University Olomouc, Šlechtitelů 27, 779 00 Olomouc, Czech Republic
| |
Collapse
|
13
|
Garcia AE. Atomistic Simulations of Thermal Unfolding. METHODS IN MOLECULAR BIOLOGY (CLIFTON, N.J.) 2022; 2376:331-341. [PMID: 34845618 DOI: 10.1007/978-1-0716-1716-8_18] [Citation(s) in RCA: 1] [Impact Index Per Article: 0.5] [Reference Citation Analysis] [Abstract] [Key Words] [Subscribe] [Scholar Register] [Indexed: 10/19/2022]
Abstract
This tutorial will provide a practical overview of the use of atomistic simulations to study thermal unfolding of biomolecules, in particular small proteins and RNA oligomers. The tutorial focuses on the use of atomistic, all atom simulations of biomolecules in explicit solvent, to study (reversible) thermal unfolding. The simulation methods described here have also been applied to study biomolecules using implicit solvent and coarse-grained models. We do not intend to provide an up-to-date review of the vast literature of biomolecular dynamics, enhanced sampling methods, force field developments, and applications of these methods. The purpose of this tutorial is to provide basic guidelines into the use of these methods to the starting scientist.
Collapse
Affiliation(s)
- Angel E Garcia
- Center for NonLinear Studies, Los Alamos National Laboratory, Los Alamos, NM, USA.
| |
Collapse
|
14
|
Ahn SH, Ojha AA, Amaro RE, McCammon JA. Gaussian-Accelerated Molecular Dynamics with the Weighted Ensemble Method: A Hybrid Method Improves Thermodynamic and Kinetic Sampling. J Chem Theory Comput 2021; 17:7938-7951. [PMID: 34844409 DOI: 10.1021/acs.jctc.1c00770] [Citation(s) in RCA: 13] [Impact Index Per Article: 4.3] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/28/2022]
Abstract
Gaussian-accelerated molecular dynamics (GaMD) is a well-established enhanced sampling method for molecular dynamics simulations that effectively samples the potential energy landscape of the system by adding a boost potential, which smoothens the surface and lowers the energy barriers between states. GaMD is unable to give time-dependent properties such as kinetics directly. On the other hand, the weighted ensemble (WE) method can efficiently sample transitions between states with its many weighted trajectories, which directly yield rates and pathways. However, convergence to equilibrium conditions remains a challenge for the WE method. Hence, we have developed a hybrid method that combines the two methods, wherein GaMD is first used to sample the potential energy landscape of the system and WE is subsequently used to further sample the potential energy landscape and kinetic properties of interest. We show that the hybrid method can sample both thermodynamic and kinetic properties more accurately and quickly compared to using either method alone.
Collapse
Affiliation(s)
- Surl-Hee Ahn
- Department of Chemistry, University of California San Diego, La Jolla 92093, California, United States
| | - Anupam A Ojha
- Department of Chemistry, University of California San Diego, La Jolla 92093, California, United States
| | - Rommie E Amaro
- Department of Chemistry, University of California San Diego, La Jolla 92093, California, United States
| | - J Andrew McCammon
- Department of Chemistry, University of California San Diego, La Jolla 92093, California, United States.,Department of Pharmacology, University of California San Diego, La Jolla 92093, California, United States
| |
Collapse
|
15
|
Unveiling the N-Terminal Homodimerization of BCL11B by Hybrid Solvent Replica-Exchange Simulations. Int J Mol Sci 2021; 22:ijms22073650. [PMID: 33807484 PMCID: PMC8036541 DOI: 10.3390/ijms22073650] [Citation(s) in RCA: 1] [Impact Index Per Article: 0.3] [Reference Citation Analysis] [Abstract] [Key Words] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 03/15/2021] [Revised: 03/29/2021] [Accepted: 03/30/2021] [Indexed: 01/28/2023] Open
Abstract
Transcription factors play a crucial role in regulating biological processes such as cell growth, differentiation, organ development and cellular signaling. Within this group, proteins equipped with zinc finger motifs (ZFs) represent the largest family of sequence-specific DNA-binding transcription regulators. Numerous studies have proven the fundamental role of BCL11B for a variety of tissues and organs such as central nervous system, T cells, skin, teeth, and mammary glands. In a previous work we identified a novel atypical zinc finger domain (CCHC-ZF) which serves as a dimerization interface of BCL11B. This domain and formation of the dimer were shown to be critically important for efficient regulation of the BCL11B target genes and could therefore represent a promising target for novel drug therapies. Here, we report the structural basis for BCL11B-BCL11B interaction mediated by the N-terminal ZF domain. By combining structure prediction algorithms, enhanced sampling molecular dynamics and fluorescence resonance energy transfer (FRET) approaches, we identified amino acid residues indispensable for the formation of the single ZF domain and directly involved in forming the dimer interface. These findings not only provide deep insight into how BCL11B acquires its active structure but also represent an important step towards rational design or selection of potential inhibitors.
Collapse
|
16
|
Zhang H, Zhang H, Chen C. Investigating the folding mechanism of the N-terminal domain of ribosomal protein L9. Proteins 2021; 89:832-844. [PMID: 33576138 DOI: 10.1002/prot.26062] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 10/07/2020] [Revised: 01/04/2021] [Accepted: 01/31/2021] [Indexed: 11/10/2022]
Abstract
Protein folding is a popular topic in the life science. However, due to the limited sampling ability of experiments and simulations, the general folding mechanism is not yet clear to us. In this work, we study the folding of the N-terminal domain of ribosomal protein L9 (NTL9) in detail by a mixing replica exchange molecular dynamics method. The simulation results are close to previous experimental observations. According to the Markov state model, the folding of the protein follows a nucleation-condensation path. Moreover, after the comparison to its 39-residue β-α-β motif, we find that the helix at the C-terminal has a great influence on the folding process of the intact protein, including the nucleation of the key residues in the transition state ensemble and the packing of the hydrophobic residues in the native state.
Collapse
Affiliation(s)
- Haozhe Zhang
- Biomolecular Physics and Modeling Group, School of Physics, Huazhong University of Science and Technology, Wuhan, China
| | - Haomiao Zhang
- Biomolecular Physics and Modeling Group, School of Physics, Huazhong University of Science and Technology, Wuhan, China
| | - Changjun Chen
- Biomolecular Physics and Modeling Group, School of Physics, Huazhong University of Science and Technology, Wuhan, China
| |
Collapse
|
17
|
Rick SW, Schwing GJ, Summa CM. An Implementation of Replica Exchange with Dynamical Scaling for Efficient Large-Scale Simulations. J Chem Inf Model 2021; 61:810-818. [PMID: 33496583 DOI: 10.1021/acs.jcim.0c01236] [Citation(s) in RCA: 4] [Impact Index Per Article: 1.3] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/30/2022]
Abstract
An implementation of the replica exchange with dynamical scaling (REDS) method in the commonly used molecular dynamics program GROMACS is presented. REDS is a replica exchange method that requires fewer replicas than conventional replica exchange while still providing data over a range of temperatures and can be used in either constant volume or constant pressure ensembles. Details for running REDS simulations are given, and an application to the human islet amyloid polypeptide (hIAPP) 11-25 fragment shows that the model efficiently samples conformational space.
Collapse
Affiliation(s)
- Steven W Rick
- Department of Chemistry, University of New Orleans, New Orleans, Louisiana 70148, United States
| | - Gregory J Schwing
- Department of Computer Science, University of New Orleans, New Orleans, Louisiana 70148, United States
| | - Christopher M Summa
- Department of Computer Science, University of New Orleans, New Orleans, Louisiana 70148, United States
| |
Collapse
|
18
|
Faizi F, Buigues PJ, Deligiannidis G, Rosta E. Simulated tempering with irreversible Gibbs sampling techniques. J Chem Phys 2020; 153:214111. [PMID: 33291930 DOI: 10.1063/5.0025775] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/14/2022] Open
Abstract
We present here two novel algorithms for simulated tempering simulations, which break the detailed balance condition (DBC) but satisfy the skewed detailed balance to ensure invariance of the target distribution. The irreversible methods we present here are based on Gibbs sampling and concern breaking DBC at the update scheme of the temperature swaps. We utilize three systems as a test bed for our methods: a Markov chain Monte Carlo simulation on a simple system described by a one-dimensional double well potential, the Ising model, and molecular dynamics simulations on alanine pentapeptide (ALA5). The relaxation times of inverse temperature, magnetic susceptibility, and energy density for the Ising model indicate clear gains in sampling efficiency over conventional Gibbs sampling techniques with DBC and also over the conventionally used simulated tempering with the Metropolis-Hastings (MH) scheme. Simulations on ALA5 with a large number of temperatures indicate distinct gains in mixing times for inverse temperature and consequently the energy of the system compared to conventional MH. With no additional computational overhead, our methods were found to be more efficient alternatives to the conventionally used simulated tempering methods with DBC. Our algorithms should be particularly advantageous in simulations of large systems with many temperature ladders, as our algorithms showed a more favorable constant scaling in Ising spin systems as compared with both reversible and irreversible MH algorithms. In future applications, our irreversible methods can also be easily tailored to utilize a given dynamical variable other than temperature to flatten rugged free energy landscapes.
Collapse
Affiliation(s)
- Fahim Faizi
- Department of Mathematics, King's College London, Strand, WC2R 2LS London, United Kingdom
| | - Pedro J Buigues
- Department of Chemistry, King's College London, 7 Trinity Street, SE1 1DB London, United Kingdom
| | - George Deligiannidis
- Department of Statistics, University of Oxford, 24-29 St Giles', OX1 3LB Oxford, United Kingdom
| | - Edina Rosta
- Department of Chemistry, King's College London, 7 Trinity Street, SE1 1DB London, United Kingdom
| |
Collapse
|
19
|
Validation of DBFOLD: An efficient algorithm for computing folding pathways of complex proteins. PLoS Comput Biol 2020; 16:e1008323. [PMID: 33196646 PMCID: PMC7704049 DOI: 10.1371/journal.pcbi.1008323] [Citation(s) in RCA: 3] [Impact Index Per Article: 0.8] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 08/26/2020] [Revised: 11/30/2020] [Accepted: 10/17/2020] [Indexed: 11/19/2022] Open
Abstract
Atomistic simulations can provide valuable, experimentally-verifiable insights into protein folding mechanisms, but existing ab initio simulation methods are restricted to only the smallest proteins due to severe computational speed limits. The folding of larger proteins has been studied using native-centric potential functions, but such models omit the potentially crucial role of non-native interactions. Here, we present an algorithm, entitled DBFOLD, which can predict folding pathways for a wide range of proteins while accounting for the effects of non-native contacts. In addition, DBFOLD can predict the relative rates of different transitions within a protein’s folding pathway. To accomplish this, rather than directly simulating folding, our method combines equilibrium Monte-Carlo simulations, which deploy enhanced sampling, with unfolding simulations at high temperatures. We show that under certain conditions, trajectories from these two types of simulations can be jointly analyzed to compute unknown folding rates from detailed balance. This requires inferring free energies from the equilibrium simulations, and extrapolating transition rates from the unfolding simulations to lower, physiologically-reasonable temperatures at which the native state is marginally stable. As a proof of principle, we show that our method can accurately predict folding pathways and Monte-Carlo rates for the well-characterized Streptococcal protein G. We then show that our method significantly reduces the amount of computation time required to compute the folding pathways of large, misfolding-prone proteins that lie beyond the reach of existing direct simulation. Our algorithm, which is available online, can generate detailed atomistic models of protein folding mechanisms while shedding light on the role of non-native intermediates which may crucially affect organismal fitness and are frequently implicated in disease. Many proteins must adopt a specific structure in order to function. Computational simulations have been used to shed light on the mechanisms of protein folding, but unfortunately, realistic simulations can typically only be run for small proteins, due to severe limits in computational speed. Here, we present a method to solve this problem, whereby instead of directly simulating folding from an unfolded state, we run simulations that allow for computation of equilibrium folding free energies, alongside high temperature simulations to compute unfolding rates. From these quantities, folding rates can be computed using detailed balance. Importantly, our method can account for the effects of nonnative contacts which transiently form during folding and must be broken prior to adoption of the native state. Such contacts, which are often excluded from simple models of folding, may crucially affect real protein folding pathways and are often observed in folding intermediates implicated in disease.
Collapse
|
20
|
Khayat E, Klimov DK, Smith AK. Phosphorylation Promotes Aβ25-35 Peptide Aggregation within the DMPC Bilayer. ACS Chem Neurosci 2020; 11:3430-3441. [PMID: 33006281 DOI: 10.1021/acschemneuro.0c00541] [Citation(s) in RCA: 8] [Impact Index Per Article: 2.0] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/30/2022] Open
Abstract
The consequences of phosphorylation of the Aβ25-35 peptide at the position Ser26 on its aggregation have not been examined. To investigate them, we performed all-atom replica exchange simulations probing the binding of phosphorylated Aβ25-35 (pAβ25-35) peptides to the dimyristoyl phosphatidylcholine (DMPC) bilayer and their subsequent aggregation. As a control, we used our previous study of unmodified peptides. We found that phosphorylation moderately reduces the helical propensity in pAβ25-35 and its binding affinity to the DMPC bilayer. Phosphorylation preserves the bimodal binding observed for unmodified Aβ25-35, which features a preferred inserted state and a less probable surface bound state. Phosphorylation also retains the inserted dimer with a head-to-tail helical aggregation interface as the most thermodynamically stable state. Importantly, this post-translation modification strengthens interpeptide interactions by adding a new aggregation "hot spot" created by cross-bridging phosphorylated Ser26 with water, cationic ions, or Lys28. The cross-bridging constitutes the molecular mechanism behind most structural phosphorylation effects. In addition, phosphorylation eliminates pAβ25-35 monomers and diversifies the pool of aggregated species. As a result, it changes the binding and aggregation mechanism by multiplying pathways leading to stable inserted dimers. These findings offer a plausible molecular rationale for experimental observations, including enhanced production of low molecular weight oligomers and cytotoxicity of phosphorylated Aβ peptides.
Collapse
Affiliation(s)
- Elias Khayat
- School of Systems Biology, George Mason University, Manassas, Virginia 20110, United States
| | - Dmitri K. Klimov
- School of Systems Biology, George Mason University, Manassas, Virginia 20110, United States
| | - Amy K. Smith
- School of Systems Biology, George Mason University, Manassas, Virginia 20110, United States
| |
Collapse
|
21
|
Kasahara K, Terazawa H, Itaya H, Goto S, Nakamura H, Takahashi T, Higo J. myPresto/omegagene 2020: a molecular dynamics simulation engine for virtual-system coupled sampling. Biophys Physicobiol 2020; 17:140-146. [PMID: 33240741 PMCID: PMC7671739 DOI: 10.2142/biophysico.bsj-2020013] [Citation(s) in RCA: 8] [Impact Index Per Article: 2.0] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 06/15/2020] [Accepted: 10/10/2020] [Indexed: 12/03/2022] Open
Abstract
The molecular dynamics (MD) method is a promising approach for investigating the molecular mechanisms of microscopic phenomena. In particular, generalized ensemble MD methods can efficiently explore the conformational space with a rugged free-energy surface. However, the implementation and acquisition of technical knowledge for each generalized ensemble MD method are not straightforward for end-users. Here, we present a new version of the myPresto/omegagene software, which is an MD simulation engine tailored for a series of generalized ensemble methods, which are virtual-system coupled multicanonical MD (V-McMD), virtual-system coupled adaptive umbrella sampling (V-AUS), and virtual-system coupled canonical MD (VcMD). This program has been applied in several studies analyzing free-energy landscapes of a variety of molecular systems with all-atom simulations. The updated version provides new functionality for coarse-grained simulations powered by the hydrophobicity scale method. The software package includes a step-by-step tutorial document for enhanced conformational sampling of the poly-glutamine (poly-Q) oligomer expressed as a one-bead per residue model. The myPresto/omegagene software is freely available at the following URL: https://github.com/kotakasahara/omegagene under the Apache2 license.
Collapse
Affiliation(s)
- Kota Kasahara
- College of Life Sciences, Ritsumeikan University, Kusatsu, Shiga 525-8577, Japan
| | - Hiroki Terazawa
- Graduate School of Life Sciences, Ritsumeikan University, Kusatsu, Shiga 525-8577, Japan
| | - Hayato Itaya
- Graduate School of Life Sciences, Ritsumeikan University, Kusatsu, Shiga 525-8577, Japan
| | - Satoshi Goto
- Graduate School of Life Sciences, Ritsumeikan University, Kusatsu, Shiga 525-8577, Japan
| | - Haruki Nakamura
- Institute for Protein Research, Osaka University, Suita, Osaka 565-0871, Japan
| | - Takuya Takahashi
- College of Life Sciences, Ritsumeikan University, Kusatsu, Shiga 525-8577, Japan
| | - Junichi Higo
- Graduate School of Simulation Studies, University of Hyogo, Kobe, Hyogo 650-0047, Japan
| |
Collapse
|
22
|
Narayan B, Yuan Y, Fathizadeh A, Elber R, Buchete NV. Long-time methods for molecular dynamics simulations: Markov State Models and Milestoning. PROGRESS IN MOLECULAR BIOLOGY AND TRANSLATIONAL SCIENCE 2020; 170:215-237. [PMID: 32145946 DOI: 10.1016/bs.pmbts.2020.01.002] [Citation(s) in RCA: 6] [Impact Index Per Article: 1.5] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Subscribe] [Scholar Register] [Indexed: 12/12/2022]
Abstract
Molecular dynamics (MD) studies of biomolecules require the ability to simulate complex biochemical systems with an increasingly larger number of particles and for longer time scales, a problem that cannot be overcome by computational hardware advances alone. A main problem springs from the intrinsically high-dimensional and complex nature of the underlying free energy landscape of most systems, and from the necessity to sample accurately such landscapes for identifying kinetic and thermodynamic states in the configurations space, and for accurate calculations of both free energy differences and of the corresponding transition rates between states. Here, we review and present applications of two increasingly popular methods that allow long-time MD simulations of biomolecular systems that can open a broad spectrum of new studies. A first approach, Markov State Models (MSMs), relies on identifying a set of configuration states in which the system resides sufficiently long to relax and loose the memory of previous transitions, and on using simulations for mapping the underlying complex energy landscape and for extracting accurate thermodynamic and kinetic information. The Markovian independence of the underlying transition probabilities creates the opportunity to increase the sampling efficiency by using sets of appropriately initialized short simulations rather than typically long MD trajectories, which also enhances sampling. This allows MSM-based studies to unveil bio-molecular mechanisms and to estimate free energy barriers with high accuracy, in a manner that is both systematic and relatively automatic, which accounts for their increasing popularity. The second approach presented, Milestoning, targets accurate studies of the ensemble of pathways connecting specific end-states (e.g., reactants and products) in a similarly systematic, accurate and highly automatic manner. Applications presented range from studies of conformational dynamics and binding of amyloid-forming peptides, cell-penetrating peptides and the DFG-flip dynamics in Abl kinase. As highlighted by the increasing number of studies using both methods, we anticipate that they will open new avenues for the investigation of systematic sampling of reactions pathways and mechanisms occurring on longer time scales than currently accessible by purely computational hardware developments.
Collapse
Affiliation(s)
- Brajesh Narayan
- School of Physics, University College Dublin, Dublin, Ireland; Institute for Discovery, University College Dublin, Dublin, Ireland
| | - Ye Yuan
- School of Physics, University College Dublin, Dublin, Ireland; Institute for Discovery, University College Dublin, Dublin, Ireland
| | - Arman Fathizadeh
- Oden Institute for Computational Engineering and Sciences, University of Texas at Austin, Austin, TX, United States
| | - Ron Elber
- Oden Institute for Computational Engineering and Sciences, University of Texas at Austin, Austin, TX, United States; Department of Chemistry, University of Texas at Austin, Austin, TX, United States
| | - Nicolae-Viorel Buchete
- School of Physics, University College Dublin, Dublin, Ireland; Institute for Discovery, University College Dublin, Dublin, Ireland.
| |
Collapse
|
23
|
Shimato T, Kasahara K, Higo J, Takahashi T. Effects of number of parallel runs and frequency of bias-strength replacement in generalized ensemble molecular dynamics simulations. PEERJ PHYSICAL CHEMISTRY 2019. [DOI: 10.7717/peerj-pchem.4] [Citation(s) in RCA: 2] [Impact Index Per Article: 0.4] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/20/2022] Open
Abstract
Background
The generalized ensemble approach with the molecular dynamics (MD) method has been widely utilized. This approach usually has two features. (i) A bias potential, whose strength is replaced during a simulation, is applied. (ii) Sampling can be performed by many parallel runs of simulations. Although the frequency of the bias-strength replacement and the number of parallel runs can be adjusted, the effects of these settings on the resultant ensemble remain unclear.
Method
In this study, we performed multicanonical MD simulations for a foldable mini-protein (Trp-cage) and two unstructured peptides (8- and 20-residue poly-glutamic acids) with various settings.
Results
As a result, running many short simulations yielded robust results for the Trp-cage model. Regarding the frequency of the bias-potential replacement, although using a high frequency enhanced the traversals in the potential energy space, it did not promote conformational changes in all the systems.
Collapse
Affiliation(s)
- Takuya Shimato
- Graduate School of Life Sciences, Ritsumeikan University, Kusatsu, Shiga, Japan
| | - Kota Kasahara
- College of Life Sciences, Ritsumeikan University, Kusatsu, Shiga, Japan
| | - Junichi Higo
- Graduate School of Simulation Studies, University of Hyogo, Kobe, Hyogo, Japan
| | - Takuya Takahashi
- College of Life Sciences, Ritsumeikan University, Kusatsu, Shiga, Japan
| |
Collapse
|
24
|
Zhang H, Gong Q, Zhang H, Chen C. Combining the biased and unbiased sampling strategy into one convenient free energy calculation method. J Comput Chem 2019; 40:1806-1815. [PMID: 30942500 DOI: 10.1002/jcc.25834] [Citation(s) in RCA: 4] [Impact Index Per Article: 0.8] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 11/28/2018] [Revised: 03/15/2019] [Accepted: 03/17/2019] [Indexed: 12/14/2022]
Abstract
Constructing a free energy landscape for a large molecule is difficult. One has to use either a high temperature or a strong driving force to enhance the sampling on the free energy barriers. In this work, we propose a mixed method that combines these two kinds of acceleration strategies into one simulation. First, it applies an adaptive biasing potential to some replicas of the molecule. These replicas are particularly accelerated in a collective variable space. Second, it places some unbiased and exchangeable replicas at various temperature levels. These replicas generate unbiased sampling data in the canonical ensemble. To improve the sampling efficiency, biased replicas transfer their state variables to the unbiased replicas after equilibrium by Monte Carlo trial moves. In comparison to previous integrated methods, it is more convenient for users. It does not need an initial reference biasing potential to guide the sampling of the molecule. And it is also unnecessary to insert many replicas for the requirement of passing the free energy barriers. The free energy calculation is accomplished in a single stage. It samples the data as fast as a biased simulation and it processes the data as simple as an unbiased simulation. The method provides a minimalist approach to the construction of the free energy landscape. © 2019 Wiley Periodicals, Inc.
Collapse
Affiliation(s)
- Haomiao Zhang
- Biomolecular Physics and Modeling Group, School of Physics, Huazhong University of Science and Technology, Wuhan 430074, Hubei, China
| | - Qiankun Gong
- Biomolecular Physics and Modeling Group, School of Physics, Huazhong University of Science and Technology, Wuhan 430074, Hubei, China
| | - Haozhe Zhang
- Biomolecular Physics and Modeling Group, School of Physics, Huazhong University of Science and Technology, Wuhan 430074, Hubei, China
| | - Changjun Chen
- Biomolecular Physics and Modeling Group, School of Physics, Huazhong University of Science and Technology, Wuhan 430074, Hubei, China
| |
Collapse
|
25
|
Kasahara K, Terazawa H, Takahashi T, Higo J. Studies on Molecular Dynamics of Intrinsically Disordered Proteins and Their Fuzzy Complexes: A Mini-Review. Comput Struct Biotechnol J 2019; 17:712-720. [PMID: 31303975 PMCID: PMC6603302 DOI: 10.1016/j.csbj.2019.06.009] [Citation(s) in RCA: 29] [Impact Index Per Article: 5.8] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 03/01/2019] [Revised: 05/29/2019] [Accepted: 06/11/2019] [Indexed: 11/19/2022] Open
Abstract
The molecular dynamics (MD) method is a promising approach toward elucidating the molecular mechanisms of intrinsically disordered regions (IDRs) of proteins and their fuzzy complexes. This mini-review introduces recent studies that apply MD simulations to investigate the molecular recognition of IDRs. Firstly, methodological issues by which MD simulations treat IDRs, such as developing force fields, treating periodic boundary conditions, and enhanced sampling approaches, are discussed. Then, several examples of the applications of MD to investigate molecular interactions of IDRs in terms of the two kinds of complex formations; coupled-folding and binding and fuzzy complex. MD simulations provide insight into the molecular mechanisms of these binding processes by sampling conformational ensembles of flexible IDRs. In particular, we focused on all-atom explicit-solvent MD simulations except for studies of higher-order assembly of IDRs. Recent advances in MD methods, and computational power make it possible to dissect the molecular details of realistic molecular systems involving the dynamic behavior of IDRs.
Collapse
Affiliation(s)
- Kota Kasahara
- College of Life Sciences, Ritsumeikan University, 1-1-1 Noji-higashi, Kusatsu, Shiga 525-8577, Japan
- Corresponding author.
| | - Hiroki Terazawa
- Graduate School of Life Sciences, Ritsumeikan University, 1-1-1 Noji-higashi, Kusatsu, Shiga 525-8577, Japan
| | - Takuya Takahashi
- College of Life Sciences, Ritsumeikan University, 1-1-1 Noji-higashi, Kusatsu, Shiga 525-8577, Japan
| | - Junichi Higo
- Graduate School of Simulation Studies, University of Hyogo, 7-1-28 Minatojima-minamimachi, Chuo-ku, Kobe 650-0047, Japan
| |
Collapse
|
26
|
Xia J, Flynn W, Gallicchio E, Uplinger K, Armstrong JD, Forli S, Olson AJ, Levy RM. Massive-Scale Binding Free Energy Simulations of HIV Integrase Complexes Using Asynchronous Replica Exchange Framework Implemented on the IBM WCG Distributed Network. J Chem Inf Model 2019; 59:1382-1397. [PMID: 30758197 PMCID: PMC6496938 DOI: 10.1021/acs.jcim.8b00817] [Citation(s) in RCA: 5] [Impact Index Per Article: 1.0] [Reference Citation Analysis] [Abstract] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/29/2022]
Abstract
To perform massive-scale replica exchange molecular dynamics (REMD) simulations for calculating binding free energies of protein-ligand complexes, we implemented the asynchronous replica exchange (AsyncRE) framework of the binding energy distribution analysis method (BEDAM) in implicit solvent on the IBM World Community Grid (WCG) and optimized the simulation parameters to reduce the overhead and improve the prediction power of the WCG AsyncRE simulations. We also performed the first massive-scale binding free energy calculations using the WCG distributed computing grid and 301 ligands from the SAMPL4 challenge for large-scale binding free energy predictions of HIV-1 integrase complexes. In total there are ∼10000 simulated complexes, ∼1 million replicas, and ∼2000 μs of aggregated MD simulations. Running AsyncRE MD simulations on the WCG requires accepting a trade-off between the number of replicas that can be run (breadth) and the number of full RE cycles that can be completed per replica (depth). As compared with synchronous Replica Exchange (SyncRE) running on tightly coupled clusters like XSEDE, on the WCG many more replicas can be launched simultaneously on heterogeneous distributed hardware, but each full RE cycle requires more overhead. We compared the WCG results with that from AutoDock and more advanced RE simulations including the use of flattening potentials to accelerate sampling of selected degrees of freedom of ligands and/or receptors related to slow dynamics due to high energy barriers. We propose a suitable strategy of RE simulations to refine high throughput docking results which can be matched to corresponding computing resources: from HPC clusters, to small or medium-size distributed campus grids, and finally to massive-scale computing networks including millions of CPUs like the resources available on the WCG.
Collapse
Affiliation(s)
- Junchao Xia
- Center for Biophysics and Computational Biology and Department of Physics , Temple University , Philadelphia , Pennsylvania 19122 , United States
| | - William Flynn
- Center for Biophysics and Computational Biology and Department of Chemistry , Temple University , Philadelphia , Pennsylvania 19122 , United States
| | - Emilio Gallicchio
- Department of Chemistry , CUNY Brooklyn College , Brooklyn , New York 11210 , United States
| | - Keith Uplinger
- IBM WCG Team, 1177 South Belt Line Road , Coppell , Texas 75019 , United States
| | - Jonathan D Armstrong
- IBM WCG Team, 11400 Burnet Road , 0453B129, Austin , Texas 78758 , United States
| | - Stefano Forli
- Department of Integrative Structural and Computational Biology , The Scripps Research Institute , La Jolla , California 92037-1000 , United States
| | - Arthur J Olson
- Department of Integrative Structural and Computational Biology , The Scripps Research Institute , La Jolla , California 92037-1000 , United States
| | - Ronald M Levy
- Center for Biophysics and Computational Biology and Department of Chemistry , Temple University , Philadelphia , Pennsylvania 19122 , United States
| |
Collapse
|
27
|
Demuynck R, Wieme J, Rogge SMJ, Dedecker KD, Vanduyfhuys L, Waroquier M, Van Speybroeck V. Protocol for Identifying Accurate Collective Variables in Enhanced Molecular Dynamics Simulations for the Description of Structural Transformations in Flexible Metal-Organic Frameworks. J Chem Theory Comput 2018; 14:5511-5526. [PMID: 30336016 PMCID: PMC6236469 DOI: 10.1021/acs.jctc.8b00725] [Citation(s) in RCA: 16] [Impact Index Per Article: 2.7] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 07/15/2018] [Indexed: 01/05/2023]
Abstract
Various kinds of flexibility have been observed in metal-organic frameworks, which may originate from the topology of the material or the presence of flexible ligands. The construction of free energy profiles describing the full dynamical behavior along the phase transition path is challenging since it is not trivial to identify collective variables able to identify all metastable states along the reaction path. In this work, a systematic three-step protocol to uniquely identify the dominant order parameters for structural transformations in flexible metal-organic frameworks and subsequently construct accurate free energy profiles is presented. Methodologically, this protocol is rooted in the time-structure based independent component analysis (tICA), a well-established statistical modeling technique embedded in the Markov state model methodology and often employed to study protein folding, that allows for the identification of the slowest order parameters characterizing the structural transformation. To ensure an unbiased and systematic identification of these order parameters, the tICA decomposition is performed based on information from a prior replica exchange (RE) simulation, as this technique enhances the sampling along all degrees of freedom of the system simultaneously. From this simulation, the tICA procedure extracts the order parameters-often structural parameters-that characterize the slowest transformations in the material. Subsequently, these order parameters are adopted in traditional enhanced sampling methods such as umbrella sampling, thermodynamic integration, and variationally enhanced sampling to construct accurate free energy profiles capturing the flexibility in these nanoporous materials. In this work, the applicability of this tICA-RE protocol is demonstrated by determining the slowest order parameters in both MIL-53(Al) and CAU-13, which exhibit a strongly different type of flexibility. The obtained free energy profiles as a function of this extracted order parameter are furthermore compared to the profiles obtained when adopting less-suited collective variables, indicating the importance of systematically selecting the relevant order parameters to construct accurate free energy profiles for flexible metal-organic frameworks, which is in correspondence with experimental findings. The method succeeds in mapping the full free energy surface in terms of appropriate collective variables for MOFs exhibiting linker flexibility. For CAU-13, we show the decreased stability of the closed pore phase by systematically adding adsorbed xylene molecules in the framework.
Collapse
Affiliation(s)
- Ruben Demuynck
- Center for Molecular Modeling, Ghent University, Technologiepark 903, B-9052 Zwijnaarde, Belgium
| | - Jelle Wieme
- Center for Molecular Modeling, Ghent University, Technologiepark 903, B-9052 Zwijnaarde, Belgium
| | - Sven M. J. Rogge
- Center for Molecular Modeling, Ghent University, Technologiepark 903, B-9052 Zwijnaarde, Belgium
| | - Karen D. Dedecker
- Center for Molecular Modeling, Ghent University, Technologiepark 903, B-9052 Zwijnaarde, Belgium
| | - Louis Vanduyfhuys
- Center for Molecular Modeling, Ghent University, Technologiepark 903, B-9052 Zwijnaarde, Belgium
| | - Michel Waroquier
- Center for Molecular Modeling, Ghent University, Technologiepark 903, B-9052 Zwijnaarde, Belgium
| | - Veronique Van Speybroeck
- Center for Molecular Modeling, Ghent University, Technologiepark 903, B-9052 Zwijnaarde, Belgium
| |
Collapse
|
28
|
Narayan B, Herbert C, Yuan Y, Rodriguez BJ, Brooks BR, Buchete NV. Conformational analysis of replica exchange MD: Temperature-dependent Markov networks for FF amyloid peptides. J Chem Phys 2018; 149:072323. [PMID: 30134732 DOI: 10.1063/1.5027580] [Citation(s) in RCA: 7] [Impact Index Per Article: 1.2] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/02/2023] Open
Abstract
Recent molecular modeling methods using Markovian descriptions of conformational states of biomolecular systems have led to powerful analysis frameworks that can accurately describe their complex dynamical behavior. In conjunction with enhanced sampling methods, such as replica exchange molecular dynamics (REMD), these frameworks allow the systematic and accurate extraction of transition probabilities between the corresponding states, in the case of Markov state models, and of statistically-optimized transition rates, in the case of the corresponding coarse master equations. However, applying automatically such methods to large molecular dynamics (MD) simulations, with explicit water molecules, remains limited both by the initial ability to identify good candidates for the underlying Markovian states and by the necessity to do so using good collective variables as reaction coordinates that allow the correct counting of inter-state transitions at various lag times. Here, we show that, in cases when representative molecular conformations can be identified for the corresponding Markovian states, and thus their corresponding collective evolution of atomic positions can be calculated along MD trajectories, one can use them to build a new type of simple collective variable, which can be particularly useful in both the correct state assignment and in the subsequent accurate counting of inter-state transition probabilities. In the case of the ubiquitously used root-mean-square deviation (RMSD) of atomic positions, we introduce the relative RMSD (RelRMSD) measure as a good reaction coordinate candidate. We apply this method to the analysis of REMD trajectories of amyloid-forming diphenylalanine (FF) peptides-a system with important nanotechnology and biomedical applications due to its self-assembling and piezoelectric properties-illustrating the use of RelRMSD in extracting its temperature-dependent intrinsic kinetics, without a priori assumptions on the functional form (e.g., Arrhenius or not) of the underlying conformational transition rates. The RelRMSD analysis enables as well a more objective assessment of the convergence of the REMD simulations. This type of collective variable may be generalized to other observables that could accurately capture conformational differences between the underlying Markov states (e.g., distance RMSD, the fraction of native contacts, etc.).
Collapse
Affiliation(s)
- Brajesh Narayan
- School of Physics, University College Dublin, Belfield, Dublin 4, Ireland
| | - Colm Herbert
- School of Physics, University College Dublin, Belfield, Dublin 4, Ireland
| | - Ye Yuan
- School of Physics, University College Dublin, Belfield, Dublin 4, Ireland
| | - Brian J Rodriguez
- School of Physics, University College Dublin, Belfield, Dublin 4, Ireland
| | - Bernard R Brooks
- Laboratory of Computational Biology, NHLBI, National Institutes of Health, Bethesda, Maryland 20892, USA
| | | |
Collapse
|
29
|
Radak BK, Suh D, Roux B. A generalized linear response framework for expanded ensemble and replica exchange simulations. J Chem Phys 2018; 149:072315. [PMID: 30134700 PMCID: PMC5984729 DOI: 10.1063/1.5027494] [Citation(s) in RCA: 7] [Impact Index Per Article: 1.2] [Reference Citation Analysis] [Abstract] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 03/02/2018] [Accepted: 05/15/2018] [Indexed: 11/15/2022] Open
Abstract
Expanded ensemble simulation is a powerful technique for enhancing sampling over a range of thermodynamic parameters. However, although the premise is relatively simple, running successful simulations in practice still presents something of an ad hoc challenge. Three main difficulties exist: (1) the selection of the thermodynamic states, (2) the selection of the sampling weights, and (3) efficient sampling of the expanded parameter space. Here we consider these problems in the context of a pairwise linear response approach to the work fluctuation theorem. The approach offers comprehensive tactics for addressing the three difficulties and can be used as either an alternative or a complement to replica exchange simulations. Importantly, the results are trivially implemented for multi-dimensional parameter spaces and they recover results from the literature aimed at the special cases of simulated/parallel tempering and replica exchange umbrella sampling. Illustrative examples are shown using the NAMD simulation engine.
Collapse
Affiliation(s)
- Brian K Radak
- Theoretical and Computational Biophysics Group, Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801-2325, USA
| | - Donghyuk Suh
- Department of Chemistry, University of Chicago, Chicago, Illinois 60637-1454, USA
| | - Benoît Roux
- Department of Biochemistry and Molecular Biology, University of Chicago, Chicago, Illinois 60637-1454, USA
| |
Collapse
|
30
|
Iwai R, Kasahara K, Takahashi T. Influence of various parameters in the replica-exchange molecular dynamics method: Number of replicas, replica-exchange frequency, and thermostat coupling time constant. Biophys Physicobiol 2018; 15:165-172. [PMID: 30250775 PMCID: PMC6145944 DOI: 10.2142/biophysico.15.0_165] [Citation(s) in RCA: 4] [Impact Index Per Article: 0.7] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 03/07/2018] [Accepted: 07/05/2018] [Indexed: 12/29/2022] Open
Abstract
The replica-exchange molecular dynamics (REMD) method has been used for conformational sampling of various biomolecular systems. To maximize sampling efficiency, some adjustable parameters must be optimized. Although it is agreed that shorter intervals between the replica-exchange attempts enhance traversals in the temperature space, details regarding the artifacts caused by these short intervals are controversial. In this study, we revisit this problem by performing REMD simulations on an alanine octapeptide in an implicit solvent. Fifty different sets of conditions, which are a combination of five replica-exchange periods, five different numbers of replicas, and two thermostat coupling time constants, were investigated. As a result, although short replica-exchange intervals enhanced the traversals in the temperature space, they led to artifacts in the ensemble average of the temperature, potential energy, and helix content. With extremely short replica-exchange intervals, i.e., attempted at every time step, the ensemble average of the temperature deviated from the thermostat temperature by ca. 7 K. Differences in the ensembles were observed even for larger replica-exchange intervals (between 100 and 1,000 steps). In addition, the shorter thermostat coupling time constant reduced the artifacts found when short replica-exchange intervals were used, implying that these artifacts are caused by insufficient thermal relaxation between the replica-exchange events. Our results will be useful to reduce the artifacts found in REMD simulations by adjusting some key parameters.
Collapse
Affiliation(s)
- Ryosuke Iwai
- Graduate School of Life Sciences, Ritsumeikan University, Kusatsu, Shiga 525-8577, Japan
| | - Kota Kasahara
- College of Life Sciences, Ritsumeikan University, Kusatsu, Shiga 525-8577, Japan
| | - Takuya Takahashi
- College of Life Sciences, Ritsumeikan University, Kusatsu, Shiga 525-8577, Japan
| |
Collapse
|
31
|
Leahy CT, Kells A, Hummer G, Buchete NV, Rosta E. Peptide dimerization-dissociation rates from replica exchange molecular dynamics. J Chem Phys 2018; 147:152725. [PMID: 29055328 DOI: 10.1063/1.5004774] [Citation(s) in RCA: 15] [Impact Index Per Article: 2.5] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/30/2022] Open
Abstract
We show how accurate rates of formation and dissociation of peptide dimers can be calculated using direct transition counting (DTC) from replica-exchange molecular dynamics (REMD) simulations. First, continuous trajectories corresponding to system replicas evolving at different temperatures are used to assign conformational states. Second, we analyze the entire REMD data to calculate the corresponding rates at each temperature directly from the number of transition counts. Finally, we compare the kinetics extracted directly, using the DTC method, with indirect estimations based on trajectory likelihood maximization using short-time propagators and on decay rates of state autocorrelation functions. For systems with relatively low-dimensional intrinsic conformational dynamics, the DTC method is simple to implement and leads to accurate temperature-dependent rates. We apply the DTC rate-extraction method to all-atom REMD simulations of dimerization of amyloid-forming NNQQ tetrapetides in explicit water. In an assessment of the REMD sampling efficiency with respect to standard MD, we find a gain of more than a factor of two at the lowest temperature.
Collapse
Affiliation(s)
- Cathal T Leahy
- School of Physics, University College Dublin, Belfield, Dublin 4, Ireland
| | - Adam Kells
- Department of Chemistry, King's College London, London SE1 1DB, United Kingdom
| | - Gerhard Hummer
- Department of Theoretical Biophysics, Max Planck Institute of Biophysics, Max-von-Laue-Straße 3, 60438 Frankfurt am Main, Germany
| | | | - Edina Rosta
- Department of Chemistry, King's College London, London SE1 1DB, United Kingdom
| |
Collapse
|
32
|
Li B, Fooksa M, Heinze S, Meiler J. Finding the needle in the haystack: towards solving the protein-folding problem computationally. Crit Rev Biochem Mol Biol 2018; 53:1-28. [PMID: 28976219 PMCID: PMC6790072 DOI: 10.1080/10409238.2017.1380596] [Citation(s) in RCA: 21] [Impact Index Per Article: 3.5] [Reference Citation Analysis] [Abstract] [Key Words] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 05/16/2017] [Revised: 08/22/2017] [Accepted: 09/13/2017] [Indexed: 12/22/2022]
Abstract
Prediction of protein tertiary structures from amino acid sequence and understanding the mechanisms of how proteins fold, collectively known as "the protein folding problem," has been a grand challenge in molecular biology for over half a century. Theories have been developed that provide us with an unprecedented understanding of protein folding mechanisms. However, computational simulation of protein folding is still difficult, and prediction of protein tertiary structure from amino acid sequence is an unsolved problem. Progress toward a satisfying solution has been slow due to challenges in sampling the vast conformational space and deriving sufficiently accurate energy functions. Nevertheless, several techniques and algorithms have been adopted to overcome these challenges, and the last two decades have seen exciting advances in enhanced sampling algorithms, computational power and tertiary structure prediction methodologies. This review aims at summarizing these computational techniques, specifically conformational sampling algorithms and energy approximations that have been frequently used to study protein-folding mechanisms or to de novo predict protein tertiary structures. We hope that this review can serve as an overview on how the protein-folding problem can be studied computationally and, in cases where experimental approaches are prohibitive, help the researcher choose the most relevant computational approach for the problem at hand. We conclude with a summary of current challenges faced and an outlook on potential future directions.
Collapse
Affiliation(s)
- Bian Li
- Department of Chemistry, Vanderbilt University, Nashville, TN, USA
- Center for Structural Biology, Vanderbilt University, Nashville, TN, USA
| | - Michaela Fooksa
- Center for Structural Biology, Vanderbilt University, Nashville, TN, USA
- Chemical and Physical Biology Graduate Program, Vanderbilt University, Nashville, TN, USA
| | - Sten Heinze
- Department of Chemistry, Vanderbilt University, Nashville, TN, USA
- Center for Structural Biology, Vanderbilt University, Nashville, TN, USA
| | - Jens Meiler
- Department of Chemistry, Vanderbilt University, Nashville, TN, USA
- Center for Structural Biology, Vanderbilt University, Nashville, TN, USA
| |
Collapse
|
33
|
Stelzl LS, Kells A, Rosta E, Hummer G. Dynamic Histogram Analysis To Determine Free Energies and Rates from Biased Simulations. J Chem Theory Comput 2017; 13:6328-6342. [PMID: 29059525 DOI: 10.1021/acs.jctc.7b00373] [Citation(s) in RCA: 47] [Impact Index Per Article: 6.7] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/28/2022]
Abstract
We present an algorithm to calculate free energies and rates from molecular simulations on biased potential energy surfaces. As input, it uses the accumulated times spent in each state or bin of a histogram and counts of transitions between them. Optimal unbiased equilibrium free energies for each of the states/bins are then obtained by maximizing the likelihood of a master equation (i.e., first-order kinetic rate model). The resulting free energies also determine the optimal rate coefficients for transitions between the states or bins on the biased potentials. Unbiased rates can be estimated, e.g., by imposing a linear free energy condition in the likelihood maximization. The resulting "dynamic histogram analysis method extended to detailed balance" (DHAMed) builds on the DHAM method. It is also closely related to the transition-based reweighting analysis method (TRAM) and the discrete TRAM (dTRAM). However, in the continuous-time formulation of DHAMed, the detailed balance constraints are more easily accounted for, resulting in compact expressions amenable to efficient numerical treatment. DHAMed produces accurate free energies in cases where the common weighted-histogram analysis method (WHAM) for umbrella sampling fails because of slow dynamics within the windows. Even in the limit of completely uncorrelated data, where WHAM is optimal in the maximum-likelihood sense, DHAMed results are nearly indistinguishable. We illustrate DHAMed with applications to ion channel conduction, RNA duplex formation, α-helix folding, and rate calculations from accelerated molecular dynamics. DHAMed can also be used to construct Markov state models from biased or replica-exchange molecular dynamics simulations. By using binless WHAM formulated as a numerical minimization problem, the bias factors for the individual states can be determined efficiently in a preprocessing step and, if needed, optimized globally afterward.
Collapse
Affiliation(s)
- Lukas S Stelzl
- Department of Theoretical Biophysics, Max Planck Institute of Biophysics , Max-von-Laue Straße 3, 60438 Frankfurt am Main, Germany
| | - Adam Kells
- Department of Chemistry, King's College London , London SE1 1UL, United Kingdom
| | - Edina Rosta
- Department of Chemistry, King's College London , London SE1 1UL, United Kingdom
| | - Gerhard Hummer
- Department of Theoretical Biophysics, Max Planck Institute of Biophysics , Max-von-Laue Straße 3, 60438 Frankfurt am Main, Germany.,Institute for Biophysics, Goethe University Frankfurt , 60438 Frankfurt am Main, Germany
| |
Collapse
|
34
|
Chen C. Constructing a multidimensional free energy surface like a spider weaving a web. J Comput Chem 2017; 38:2298-2306. [PMID: 28718973 DOI: 10.1002/jcc.24881] [Citation(s) in RCA: 1] [Impact Index Per Article: 0.1] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 03/31/2017] [Revised: 06/20/2017] [Accepted: 06/22/2017] [Indexed: 01/13/2023]
Abstract
Complete free energy surface in the collective variable space provides important information of the reaction mechanisms of the molecules. But, sufficient sampling in the collective variable space is not easy. The space expands quickly with the number of the collective variables. To solve the problem, many methods utilize artificial biasing potentials to flatten out the original free energy surface of the molecule in the simulation. Their performances are sensitive to the definitions of the biasing potentials. Fast-growing biasing potential accelerates the sampling speed but decreases the accuracy of the free energy result. Slow-growing biasing potential gives an optimized result but needs more simulation time. In this article, we propose an alternative method. It adds the biasing potential to a representative point of the molecule in the collective variable space to improve the conformational sampling. And the free energy surface is calculated from the free energy gradient in the constrained simulation, not given by the negative of the biasing potential as previous methods. So the presented method does not require the biasing potential to remove all the barriers and basins on the free energy surface exactly. Practical applications show that the method in this work is able to produce the accurate free energy surfaces for different molecules in a short time period. The free energy errors are small in the cases of various biasing potentials. © 2017 Wiley Periodicals, Inc.
Collapse
Affiliation(s)
- Changjun Chen
- Biomolecular Physics and Modeling Group, School of Physics, Huazhong University of Science and Technology, Wuhan, Hubei, 430074, China
| |
Collapse
|
35
|
Stelzl LS, Hummer G. Kinetics from Replica Exchange Molecular Dynamics Simulations. J Chem Theory Comput 2017; 13:3927-3935. [DOI: 10.1021/acs.jctc.7b00372] [Citation(s) in RCA: 50] [Impact Index Per Article: 7.1] [Reference Citation Analysis] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/28/2022]
Affiliation(s)
- Lukas S. Stelzl
- Department
of Theoretical Biophysics, Max Planck Institute of Biophysics, Max-von-Laue
Straße 3, 60438 Frankfurt am Main, Germany
| | - Gerhard Hummer
- Department
of Theoretical Biophysics, Max Planck Institute of Biophysics, Max-von-Laue
Straße 3, 60438 Frankfurt am Main, Germany
- Institute
for Biophysics, Goethe University Frankfurt, 60438 Frankfurt
am Main, Germany
| |
Collapse
|
36
|
Islam B, Stadlbauer P, Gil-Ley A, Pérez-Hernández G, Haider S, Neidle S, Bussi G, Banas P, Otyepka M, Sponer J. Exploring the Dynamics of Propeller Loops in Human Telomeric DNA Quadruplexes Using Atomistic Simulations. J Chem Theory Comput 2017; 13:2458-2480. [PMID: 28475322 PMCID: PMC5514396 DOI: 10.1021/acs.jctc.7b00226] [Citation(s) in RCA: 33] [Impact Index Per Article: 4.7] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/24/2023]
Abstract
![]()
We
have carried out a series of extended unbiased molecular dynamics
(MD) simulations (up to 10 μs long, ∼162 μs in
total) complemented by replica-exchange with the collective variable
tempering (RECT) approach for several human telomeric DNA G-quadruplex
(GQ) topologies with TTA propeller loops. We used different AMBER
DNA force-field variants and also processed simulations by Markov
State Model (MSM) analysis. The slow conformational transitions in
the propeller loops took place on a scale of a few μs, emphasizing
the need for long simulations in studies of GQ dynamics. The propeller
loops sampled similar ensembles for all GQ topologies and for all
force-field dihedral-potential variants. The outcomes of standard
and RECT simulations were consistent and captured similar spectrum
of loop conformations. However, the most common crystallographic loop
conformation was very unstable with all force-field versions. Although
the loss of canonical γ-trans state of the
first propeller loop nucleotide could be related to the indispensable
bsc0 α/γ dihedral potential, even supporting this particular
dihedral by a bias was insufficient to populate the experimentally
dominant loop conformation. In conclusion, while our simulations were
capable of providing a reasonable albeit not converged sampling of
the TTA propeller loop conformational space, the force-field description
still remained far from satisfactory.
Collapse
Affiliation(s)
- Barira Islam
- Institute of Biophysics, Academy of Sciences of the Czech Republic , Královopolská 135, 612 65 Brno, Czech Republic
| | - Petr Stadlbauer
- Institute of Biophysics, Academy of Sciences of the Czech Republic , Královopolská 135, 612 65 Brno, Czech Republic.,Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Faculty of Science, Palacký University , 17. listopadu 1192/12, 771 46 Olomouc, Czech Republic
| | - Alejandro Gil-Ley
- Scuola Internazionale Superiore di Studi Avanzati, Via Bonomea 265, 34136 Trieste, Italy
| | - Guillermo Pérez-Hernández
- Department for Mathematics and Computer Science, Freie Universität Berlin , Arnimallee 6, Berlin 14195, Germany
| | - Shozeb Haider
- UCL School of Pharmacy, 29-39 Brunswick Square, London WC1N 1AX, U.K
| | - Stephen Neidle
- UCL School of Pharmacy, 29-39 Brunswick Square, London WC1N 1AX, U.K
| | - Giovanni Bussi
- Scuola Internazionale Superiore di Studi Avanzati, Via Bonomea 265, 34136 Trieste, Italy
| | - Pavel Banas
- Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Faculty of Science, Palacký University , 17. listopadu 1192/12, 771 46 Olomouc, Czech Republic
| | - Michal Otyepka
- Institute of Biophysics, Academy of Sciences of the Czech Republic , Královopolská 135, 612 65 Brno, Czech Republic.,Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Faculty of Science, Palacký University , 17. listopadu 1192/12, 771 46 Olomouc, Czech Republic
| | - Jiri Sponer
- Institute of Biophysics, Academy of Sciences of the Czech Republic , Královopolská 135, 612 65 Brno, Czech Republic.,Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Faculty of Science, Palacký University , 17. listopadu 1192/12, 771 46 Olomouc, Czech Republic
| |
Collapse
|
37
|
Tan Z, Xia J, Zhang BW, Levy RM. Locally weighted histogram analysis and stochastic solution for large-scale multi-state free energy estimation. J Chem Phys 2016; 144:034107. [PMID: 26801020 DOI: 10.1063/1.4939768] [Citation(s) in RCA: 14] [Impact Index Per Article: 1.8] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/28/2022] Open
Abstract
The weighted histogram analysis method (WHAM) including its binless extension has been developed independently in several different contexts, and widely used in chemistry, physics, and statistics, for computing free energies and expectations from multiple ensembles. However, this method, while statistically efficient, is computationally costly or even infeasible when a large number, hundreds or more, of distributions are studied. We develop a locally WHAM (local WHAM) from the perspective of simulations of simulations (SOS), using generalized serial tempering (GST) to resample simulated data from multiple ensembles. The local WHAM equations based on one jump attempt per GST cycle can be solved by optimization algorithms orders of magnitude faster than standard implementations of global WHAM, but yield similarly accurate estimates of free energies to global WHAM estimates. Moreover, we propose an adaptive SOS procedure for solving local WHAM equations stochastically when multiple jump attempts are performed per GST cycle. Such a stochastic procedure can lead to more accurate estimates of equilibrium distributions than local WHAM with one jump attempt per cycle. The proposed methods are broadly applicable when the original data to be "WHAMMED" are obtained properly by any sampling algorithm including serial tempering and parallel tempering (replica exchange). To illustrate the methods, we estimated absolute binding free energies and binding energy distributions using the binding energy distribution analysis method from one and two dimensional replica exchange molecular dynamics simulations for the beta-cyclodextrin-heptanoate host-guest system. In addition to the computational advantage of handling large datasets, our two dimensional WHAM analysis also demonstrates that accurate results similar to those from well-converged data can be obtained from simulations for which sampling is limited and not fully equilibrated.
Collapse
Affiliation(s)
- Zhiqiang Tan
- Department of Statistics, Rutgers University, Piscataway, New Jersey 08854, USA
| | - Junchao Xia
- Center for Biophysics and Computational Biology, Department of Chemistry and Institute for Computational Molecular Science, Temple University, Philadelphia, Pennsylvania 19122, USA
| | - Bin W Zhang
- Center for Biophysics and Computational Biology, Department of Chemistry and Institute for Computational Molecular Science, Temple University, Philadelphia, Pennsylvania 19122, USA
| | - Ronald M Levy
- Center for Biophysics and Computational Biology, Department of Chemistry and Institute for Computational Molecular Science, Temple University, Philadelphia, Pennsylvania 19122, USA
| |
Collapse
|
38
|
Yu TQ, Lu J, Abrams CF, Vanden-Eijnden E. Multiscale implementation of infinite-swap replica exchange molecular dynamics. Proc Natl Acad Sci U S A 2016; 113:11744-11749. [PMID: 27698148 PMCID: PMC5081654 DOI: 10.1073/pnas.1605089113] [Citation(s) in RCA: 18] [Impact Index Per Article: 2.3] [Reference Citation Analysis] [Abstract] [Key Words] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/18/2022] Open
Abstract
Replica exchange molecular dynamics (REMD) is a popular method to accelerate conformational sampling of complex molecular systems. The idea is to run several replicas of the system in parallel at different temperatures that are swapped periodically. These swaps are typically attempted every few MD steps and accepted or rejected according to a Metropolis-Hastings criterion. This guarantees that the joint distribution of the composite system of replicas is the normalized sum of the symmetrized product of the canonical distributions of these replicas at the different temperatures. Here we propose a different implementation of REMD in which (i) the swaps obey a continuous-time Markov jump process implemented via Gillespie's stochastic simulation algorithm (SSA), which also samples exactly the aforementioned joint distribution and has the advantage of being rejection free, and (ii) this REMD-SSA is combined with the heterogeneous multiscale method to accelerate the rate of the swaps and reach the so-called infinite-swap limit that is known to optimize sampling efficiency. The method is easy to implement and can be trivially parallelized. Here we illustrate its accuracy and efficiency on the examples of alanine dipeptide in vacuum and C-terminal β-hairpin of protein G in explicit solvent. In this latter example, our results indicate that the landscape of the protein is a triple funnel with two folded structures and one misfolded structure that are stabilized by H-bonds.
Collapse
Affiliation(s)
- Tang-Qing Yu
- Courant Institute of Mathematical Sciences, New York University, New York, NY 10012
| | - Jianfeng Lu
- Department of Mathematics, Duke University, Durham, NC 27708; Department of Physics, Duke University, Durham, NC 27708; Department of Chemistry, Duke University, Durham, NC 27708
| | - Cameron F Abrams
- Department of Chemical and Biological Engineering, Drexel University, Philadelphia, PA 19104
| | - Eric Vanden-Eijnden
- Courant Institute of Mathematical Sciences, New York University, New York, NY 10012;
| |
Collapse
|
39
|
Pickard FC, König G, Simmonett AC, Shao Y, Brooks BR. An efficient protocol for obtaining accurate hydration free energies using quantum chemistry and reweighting from molecular dynamics simulations. Bioorg Med Chem 2016; 24:4988-4997. [PMID: 27667551 DOI: 10.1016/j.bmc.2016.08.031] [Citation(s) in RCA: 13] [Impact Index Per Article: 1.6] [Reference Citation Analysis] [Abstract] [Key Words] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 06/16/2016] [Revised: 08/18/2016] [Accepted: 08/19/2016] [Indexed: 11/15/2022]
Abstract
The non-Boltzmann Bennett (NBB) free energy estimator method is applied to 21 molecules from the blind subset of the SAMPL4 challenge. When NBB is applied with the SMD implicit solvent model, and the OLYP/DZP level of quantum chemistry, highly accurate hydration free energy calculations are obtained with respect to experiment (RMSD=0.89kcal·mol-1). Other quantum chemical methods are also tested, and the effects of solvent model, density functional, basis set are explored in this benchmarking study, providing a framework for improvements in calculating hydration free energies. We provide a practical guide for using the best QM-NBB protocols that are consistently more accurate than either pure QM or pure MM alone. In situations where high accuracy hydration free energy predictions are needed, the QM-NBB method with SMD implicit solvent should be the first choice of quantum chemists.
Collapse
Affiliation(s)
- Frank C Pickard
- National Institutes of Health - National Heart, Lung and Blood Institute, Laboratory of Computational Biology, 5635 Fishers Lane, T-900 Suite, Rockville, MD 20852, USA
| | - Gerhard König
- National Institutes of Health - National Heart, Lung and Blood Institute, Laboratory of Computational Biology, 5635 Fishers Lane, T-900 Suite, Rockville, MD 20852, USA; Max Planck Institut für Kohlenforschung, 45470 Mülheim an der Ruhr, NRW, Germany
| | - Andrew C Simmonett
- National Institutes of Health - National Heart, Lung and Blood Institute, Laboratory of Computational Biology, 5635 Fishers Lane, T-900 Suite, Rockville, MD 20852, USA
| | - Yihan Shao
- Department of Chemistry and Biochemistry University of Oklahoma Norman, OK 73019, USA
| | - Bernard R Brooks
- National Institutes of Health - National Heart, Lung and Blood Institute, Laboratory of Computational Biology, 5635 Fishers Lane, T-900 Suite, Rockville, MD 20852, USA
| |
Collapse
|
40
|
Leahy CT, Murphy RD, Hummer G, Rosta E, Buchete NV. Coarse Master Equations for Binding Kinetics of Amyloid Peptide Dimers. J Phys Chem Lett 2016; 7:2676-2682. [PMID: 27323250 DOI: 10.1021/acs.jpclett.6b00518] [Citation(s) in RCA: 24] [Impact Index Per Article: 3.0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 06/06/2023]
Abstract
We characterize the kinetics of dimer formation of the short amyloid microcrystal-forming tetrapeptides NNQQ by constructing coarse master equations for the conformational dynamics of the system, using temperature replica-exchange molecular dynamics (REMD) simulations. We minimize the effects of Kramers-type recrossings by assigning conformational states based on their sequential time evolution. Transition rates are further estimated from short-time state propagators by maximizing the likelihood that the extracted rates agree with the observed atomistic trajectories without any a priori assumptions about their temperature dependence. Here, we evaluate the rates for both continuous replica trajectories that visit different temperatures and for discontinuous data corresponding to each REMD temperature. While the binding-unbinding kinetic process is clearly Markovian, the conformational dynamics of the bound NNQQ dimer has a complex character. Our kinetic analysis allows us to discriminate between short-lived encounter pairs and strongly bound conformational states. The conformational dynamics of NNQQ dimers supports a kinetically driven aggregation mechanism, in agreement with the polymorphic character reported for amyloid aggregates such as microcrystals and fibrils.
Collapse
Affiliation(s)
- Cathal T Leahy
- School of Physics, University College Dublin , Belfield, Dublin 4, Ireland
- Complex and Adaptive Systems Laboratory, University College Dublin , Belfield, Dublin 4, Ireland
| | - Ronan D Murphy
- School of Physics, University College Dublin , Belfield, Dublin 4, Ireland
- Complex and Adaptive Systems Laboratory, University College Dublin , Belfield, Dublin 4, Ireland
| | - Gerhard Hummer
- Department of Theoretical Biophysics, Max Planck Institute of Biophysics , Max-von-Laue-Straße 3, D-60438 Frankfurt am Main, Germany
| | - Edina Rosta
- Department of Chemistry, King's College London , London SE1 1DB, United Kingdom
| | - Nicolae-Viorel Buchete
- School of Physics, University College Dublin , Belfield, Dublin 4, Ireland
- Complex and Adaptive Systems Laboratory, University College Dublin , Belfield, Dublin 4, Ireland
| |
Collapse
|
41
|
Zhang BW, Dai W, Gallicchio E, He P, Xia J, Tan Z, Levy RM. Simulating Replica Exchange: Markov State Models, Proposal Schemes, and the Infinite Swapping Limit. J Phys Chem B 2016; 120:8289-301. [PMID: 27079355 DOI: 10.1021/acs.jpcb.6b02015] [Citation(s) in RCA: 29] [Impact Index Per Article: 3.6] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/18/2023]
Abstract
Replica exchange molecular dynamics is a multicanonical simulation technique commonly used to enhance the sampling of solvated biomolecules on rugged free energy landscapes. While replica exchange is relatively easy to implement, there are many unanswered questions about how to use this technique most efficiently, especially because it is frequently the case in practice that replica exchange simulations are not fully converged. A replica exchange cycle consists of a series of molecular dynamics steps of a set of replicas moving under different Hamiltonians or at different thermodynamic states followed by one or more replica exchange attempts to swap replicas among the different states. How the replica exchange cycle is constructed affects how rapidly the system equilibrates. We have constructed a Markov state model of replica exchange (MSMRE) using long molecular dynamics simulations of a host-guest binding system as an example, in order to study how different implementations of the replica exchange cycle can affect the sampling efficiency. We analyze how the number of replica exchange attempts per cycle, the number of MD steps per cycle, and the interaction between the two parameters affects the largest implied time scale of the MSMRE simulation. The infinite swapping limit is an important concept in replica exchange. We show how to estimate the infinite swapping limit from the diagonal elements of the exchange transition matrix constructed from MSMRE "simulations of simulations" as well as from relatively short runs of the actual replica exchange simulations.
Collapse
Affiliation(s)
- Bin W Zhang
- Center for Biophysics and Computational Biology, Department of Chemistry and Institute for Computational Molecular Science, Temple University , Philadelphia, Pennsylvania 19122, United States
| | - Wei Dai
- Department of Physics and Astronomy, Rutgers, the State University of New Jersey , Piscataway, New Jersey 08854, United States
| | - Emilio Gallicchio
- Department of Chemistry, Brooklyn College of the City University of New York , Brooklyn, New York 11210, United States
| | - Peng He
- Center for Biophysics and Computational Biology, Department of Chemistry and Institute for Computational Molecular Science, Temple University , Philadelphia, Pennsylvania 19122, United States
| | - Junchao Xia
- Center for Biophysics and Computational Biology, Department of Chemistry and Institute for Computational Molecular Science, Temple University , Philadelphia, Pennsylvania 19122, United States
| | - Zhiqiang Tan
- Department of Statistics, Rutgers, the State University of New Jersey , Piscataway, New Jersey 08854, United States
| | - Ronald M Levy
- Center for Biophysics and Computational Biology, Department of Chemistry and Institute for Computational Molecular Science, Temple University , Philadelphia, Pennsylvania 19122, United States
| |
Collapse
|
42
|
Chen C, Huang Y. Walking freely in the energy and temperature space by the modified replica exchange molecular dynamics method. J Comput Chem 2016; 37:1565-75. [DOI: 10.1002/jcc.24371] [Citation(s) in RCA: 4] [Impact Index Per Article: 0.5] [Reference Citation Analysis] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 01/10/2016] [Revised: 03/01/2016] [Accepted: 03/02/2016] [Indexed: 11/05/2022]
Affiliation(s)
- Changjun Chen
- Biomolecular Physics and Modelling Group, School of Physics; Huazhong University of Science and Technology; Wuhan Hubei 430074 China
| | - Yanzhao Huang
- Biomolecular Physics and Modelling Group, School of Physics; Huazhong University of Science and Technology; Wuhan Hubei 430074 China
| |
Collapse
|
43
|
Neale C, Pomès R, García AE. Peptide Bond Isomerization in High-Temperature Simulations. J Chem Theory Comput 2016; 12:1989-99. [PMID: 26866899 DOI: 10.1021/acs.jctc.5b01022] [Citation(s) in RCA: 17] [Impact Index Per Article: 2.1] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/29/2022]
Abstract
Force fields for molecular simulation are generally optimized to model macromolecules such as proteins at ambient temperature and pressure. Nevertheless, elevated temperatures are frequently used to enhance conformational sampling, either during system setup or as a component of an advanced sampling technique such as temperature replica exchange. Because macromolecular force fields are now put upon to simulate temperatures and time scales that greatly exceed their original design specifications, it is appropriate to re-evaluate whether these force fields are up to the task. Here, we quantify the rates of peptide bond isomerization in high-temperature simulations of three octameric peptides and a small fast-folding protein. We show that peptide octamers with and without proline residues undergo cis/trans isomerization every 1-5 ns at 800 K with three classical atomistic force fields (AMBER99SB-ILDN, CHARMM22/CMAP, and OPLS-AA/L). On the low microsecond time scale, these force fields permit isomerization of nonprolyl peptide bonds at temperatures ≥500 K, and the CHARMM22/CMAP force field permits isomerization of prolyl peptide bonds ≥400 K. Moreover, the OPLS-AA/L force field allows chiral inversion about the Cα atom at 800 K. Finally, we show that temperature replica exchange permits cis peptide bonds developed at 540 K to subsequently migrate back to the 300 K ensemble, where cis peptide bonds are present in 2 ± 1% of the population of Trp-cage TC5b, including up to 4% of its folded state. Further work is required to assess the accuracy of cis/trans isomerization in the current generation of protein force fields.
Collapse
Affiliation(s)
- Chris Neale
- Center for NonLinear Studies (CNLS), MS B258, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, United States
| | - Régis Pomès
- Molecular Structure and Function, The Hospital for Sick Children , 686 Bay Street, Toronto, Ontario M5G 0A4, Canada.,Department of Biochemistry, University of Toronto , 101 College Street, Toronto, Ontario M5G 1L7, Canada
| | - Angel E García
- Center for NonLinear Studies (CNLS), MS B258, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, United States
| |
Collapse
|
44
|
Pan AC, Weinreich TM, Piana S, Shaw DE. Demonstrating an Order-of-Magnitude Sampling Enhancement in Molecular Dynamics Simulations of Complex Protein Systems. J Chem Theory Comput 2016; 12:1360-7. [PMID: 26866996 DOI: 10.1021/acs.jctc.5b00913] [Citation(s) in RCA: 61] [Impact Index Per Article: 7.6] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/29/2022]
Abstract
Molecular dynamics (MD) simulations can describe protein motions in atomic detail, but transitions between protein conformational states sometimes take place on time scales that are infeasible or very expensive to reach by direct simulation. Enhanced sampling methods, the aim of which is to increase the sampling efficiency of MD simulations, have thus been extensively employed. The effectiveness of such methods when applied to complex biological systems like proteins, however, has been difficult to establish because even enhanced sampling simulations of such systems do not typically reach time scales at which convergence is extensive enough to reliably quantify sampling efficiency. Here, we obtain sufficiently converged simulations of three proteins to evaluate the performance of simulated tempering, a member of a widely used class of enhanced sampling methods that use elevated temperature to accelerate sampling. Simulated tempering simulations with individual lengths of up to 100 μs were compared to (previously published) conventional MD simulations with individual lengths of up to 1 ms. With two proteins, BPTI and ubiquitin, we evaluated the efficiency of sampling of conformational states near the native state, and for the third, the villin headpiece, we examined the rate of folding and unfolding. Our comparisons demonstrate that simulated tempering can consistently achieve a substantial sampling speedup of an order of magnitude or more relative to conventional MD.
Collapse
Affiliation(s)
- Albert C Pan
- D. E. Shaw Research , New York, New York 10036, United States
| | | | - Stefano Piana
- D. E. Shaw Research , New York, New York 10036, United States
| | - David E Shaw
- D. E. Shaw Research , New York, New York 10036, United States.,Department of Biochemistry and Molecular Biophysics, Columbia University , New York, New York 10032, United States
| |
Collapse
|
45
|
The good, the bad and the user in soft matter simulations. BIOCHIMICA ET BIOPHYSICA ACTA-BIOMEMBRANES 2016; 1858:2529-2538. [PMID: 26862882 DOI: 10.1016/j.bbamem.2016.02.004] [Citation(s) in RCA: 74] [Impact Index Per Article: 9.3] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Subscribe] [Scholar Register] [Received: 12/16/2015] [Revised: 02/01/2016] [Accepted: 02/02/2016] [Indexed: 11/21/2022]
Abstract
Molecular dynamics (MD) simulations have become popular in materials science, biochemistry, biophysics and several other fields. Improvements in computational resources, in quality of force field parameters and algorithms have yielded significant improvements in performance and reliability. On the other hand, no method of research is error free. In this review, we discuss a few examples of errors and artifacts due to various sources and discuss how to avoid them. Besides bringing attention to artifacts and proper practices in simulations, we also aim to provide the reader with a starting point to explore these issues further. In particular, we hope that the discussion encourages researchers to check software, parameters, protocols and, most importantly, their own practices in order to minimize the possibility of errors. The focus here is on practical issues. This article is part of a Special Issue entitled: Biosimulations edited by Ilpo Vattulainen and Tomasz Róg.
Collapse
|
46
|
Enhanced, targeted sampling of high-dimensional free-energy landscapes using variationally enhanced sampling, with an application to chignolin. Proc Natl Acad Sci U S A 2016; 113:1150-5. [PMID: 26787868 DOI: 10.1073/pnas.1519712113] [Citation(s) in RCA: 41] [Impact Index Per Article: 5.1] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/18/2022] Open
Abstract
The capabilities of molecular simulations have been greatly extended by a number of widely used enhanced sampling methods that facilitate escaping from metastable states and crossing large barriers. Despite these developments there are still many problems which remain out of reach for these methods which has led to a vigorous effort in this area. One of the most important problems that remains unsolved is sampling high-dimensional free-energy landscapes and systems that are not easily described by a small number of collective variables. In this work we demonstrate a new way to compute free-energy landscapes of high dimensionality based on the previously introduced variationally enhanced sampling, and we apply it to the miniprotein chignolin.
Collapse
|
47
|
Xia J, Flynn WF, Gallicchio E, Zhang BW, He P, Tan Z, Levy RM. Large-scale asynchronous and distributed multidimensional replica exchange molecular simulations and efficiency analysis. J Comput Chem 2015; 36:1772-85. [PMID: 26149645 PMCID: PMC4512903 DOI: 10.1002/jcc.23996] [Citation(s) in RCA: 18] [Impact Index Per Article: 2.0] [Reference Citation Analysis] [Abstract] [Key Words] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 02/20/2015] [Revised: 06/10/2015] [Accepted: 06/11/2015] [Indexed: 01/25/2023]
Abstract
We describe methods to perform replica exchange molecular dynamics (REMD) simulations asynchronously (ASyncRE). The methods are designed to facilitate large scale REMD simulations on grid computing networks consisting of heterogeneous and distributed computing environments as well as on homogeneous high-performance clusters. We have implemented these methods on NSF (National Science Foundation) XSEDE (Extreme Science and Engineering Discovery Environment) clusters and BOINC (Berkeley Open Infrastructure for Network Computing) distributed computing networks at Temple University and Brooklyn College at CUNY (the City University of New York). They are also being implemented on the IBM World Community Grid. To illustrate the methods, we have performed extensive (more than 60 ms in aggregate) simulations for the beta-cyclodextrin-heptanoate host-guest system in the context of one- and two-dimensional ASyncRE, and we used the results to estimate absolute binding free energies using the binding energy distribution analysis method. We propose ways to improve the efficiency of REMD simulations: these include increasing the number of exchanges attempted after a specified molecular dynamics (MD) period up to the fast exchange limit and/or adjusting the MD period to allow sufficient internal relaxation within each thermodynamic state. Although ASyncRE simulations generally require long MD periods (>picoseconds) per replica exchange cycle to minimize the overhead imposed by heterogeneous computing networks, we found that it is possible to reach an efficiency similar to conventional synchronous REMD, by optimizing the combination of the MD period and the number of exchanges attempted per cycle.
Collapse
Affiliation(s)
- Junchao Xia
- Center for Biophysics and Computational Biology, Department of Chemistry and Institute for Computational Molecular Science, Temple University, Philadelphia, PA 19122
| | - William F. Flynn
- Center for Biophysics and Computational Biology, Department of Chemistry and Institute for Computational Molecular Science, Temple University, Philadelphia, PA 19122
- Department of Physics & Astronomy, Rutgers University, Piscataway, NJ 08854
| | | | - Bin W. Zhang
- Center for Biophysics and Computational Biology, Department of Chemistry and Institute for Computational Molecular Science, Temple University, Philadelphia, PA 19122
| | - Peng He
- Center for Biophysics and Computational Biology, Department of Chemistry and Institute for Computational Molecular Science, Temple University, Philadelphia, PA 19122
| | - Zhiqiang Tan
- Department of Statistics, Rutgers University, Piscataway, NJ 08854
| | - Ronald M. Levy
- Center for Biophysics and Computational Biology, Department of Chemistry and Institute for Computational Molecular Science, Temple University, Philadelphia, PA 19122
| |
Collapse
|
48
|
Chen C, Xiao Y, Huang Y. Improving the replica-exchange molecular-dynamics method for efficient sampling in the temperature space. PHYSICAL REVIEW. E, STATISTICAL, NONLINEAR, AND SOFT MATTER PHYSICS 2015; 91:052708. [PMID: 26066200 DOI: 10.1103/physreve.91.052708] [Citation(s) in RCA: 3] [Impact Index Per Article: 0.3] [Reference Citation Analysis] [Abstract] [MESH Headings] [Track Full Text] [Subscribe] [Scholar Register] [Received: 01/06/2015] [Indexed: 06/04/2023]
Abstract
Replica-exchange molecular dynamics (REMD) is a popular sampling method in the molecular simulation. By frequently exchanging the replicas at different temperatures, the molecule can jump out of the minima and sample efficiently in the conformational space. Although REMD has been shown to be practical in a lot of applications, it does have a critical limitation. All the replicas at all the temperatures must be simulated for a period between the replica-exchange steps. This may be problematic for the reaction with high free energy barriers. In that case, too many replicas are required in the simulation. To reduce the calculation quantity and improve its performance, in this paper we propose a modified REMD method. During the simulation, each replica at each temperature can stay in either the active or inactive state and only switch between the states at the exchange step. In the active state, the replica moves freely in the canonical ensemble by the normal molecular dynamics, and in the inactive state, the replica is frozen temporarily until the next exchange step. The number of the replicas in the active states (active replicas) depends on the number of CPUs in the computer. Using the additional inactive replicas, one can perform an REMD simulation in a wider temperature space. The practical applications show that the modified REMD method is reliable. With the same number of active replicas, this REMD method can produce a more reasonable free energy surface around the free energy minima than the standard REMD method.
Collapse
Affiliation(s)
- Changjun Chen
- Biomolecular Physics and Modelling Group, School of Physics, Huazhong University of Science and Technology, Wuhan 430074, Hubei, China
| | - Yi Xiao
- Biomolecular Physics and Modelling Group, School of Physics, Huazhong University of Science and Technology, Wuhan 430074, Hubei, China
| | - Yanzhao Huang
- Biomolecular Physics and Modelling Group, School of Physics, Huazhong University of Science and Technology, Wuhan 430074, Hubei, China
| |
Collapse
|
49
|
Sharir-Ivry A, Varatharaj R, Shurki A. Challenges within the Linear Response Approximation When Studying Enzyme Catalysis and Effects of Mutations. J Chem Theory Comput 2014; 11:293-302. [DOI: 10.1021/ct500751f] [Citation(s) in RCA: 6] [Impact Index Per Article: 0.6] [Reference Citation Analysis] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/19/2022]
Affiliation(s)
- Avital Sharir-Ivry
- Institute for Drug Design,
School of Pharmacy, The Lise Meitner-Minerva Center for Computational
Quantum Chemistry, The Hebrew University of Jerusalem, Jerusalem 91120, Israel
| | - Rajapandian Varatharaj
- Institute for Drug Design,
School of Pharmacy, The Lise Meitner-Minerva Center for Computational
Quantum Chemistry, The Hebrew University of Jerusalem, Jerusalem 91120, Israel
| | - Avital Shurki
- Institute for Drug Design,
School of Pharmacy, The Lise Meitner-Minerva Center for Computational
Quantum Chemistry, The Hebrew University of Jerusalem, Jerusalem 91120, Israel
| |
Collapse
|
50
|
Bacci M, Vitalis A, Caflisch A. A molecular simulation protocol to avoid sampling redundancy and discover new states. Biochim Biophys Acta Gen Subj 2014; 1850:889-902. [PMID: 25193737 DOI: 10.1016/j.bbagen.2014.08.013] [Citation(s) in RCA: 18] [Impact Index Per Article: 1.8] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 06/18/2014] [Revised: 08/21/2014] [Accepted: 08/25/2014] [Indexed: 10/24/2022]
Abstract
BACKGROUND For biomacromolecules or their assemblies, experimental knowledge is often restricted to specific states. Ambiguity pervades simulations of these complex systems because there is no prior knowledge of relevant phase space domains, and sampling recurrence is difficult to achieve. In molecular dynamics methods, ruggedness of the free energy surface exacerbates this problem by slowing down the unbiased exploration of phase space. Sampling is inefficient if dwell times in metastable states are large. METHODS We suggest a heuristic algorithm to terminate and reseed trajectories run in multiple copies in parallel. It uses a recent method to order snapshots, which provides notions of "interesting" and "unique" for individual simulations. We define criteria to guide the reseeding of runs from more "interesting" points if they sample overlapping regions of phase space. RESULTS Using a pedagogical example and an α-helical peptide, the approach is demonstrated to amplify the rate of exploration of phase space and to discover metastable states not found by conventional sampling schemes. Evidence is provided that accurate kinetics and pathways can be extracted from the simulations. CONCLUSIONS The method, termed PIGS for Progress Index Guided Sampling, proceeds in unsupervised fashion, is scalable, and benefits synergistically from larger numbers of replicas. Results confirm that the underlying ideas are appropriate and sufficient to enhance sampling. GENERAL SIGNIFICANCE In molecular simulations, errors caused by not exploring relevant domains in phase space are always unquantifiable and can be arbitrarily large. Our protocol adds to the toolkit available to researchers in reducing these types of errors. This article is part of a Special Issue entitled "Recent developments of molecular dynamics".
Collapse
Affiliation(s)
- Marco Bacci
- University of Zurich, Department of Biochemistry, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland
| | - Andreas Vitalis
- University of Zurich, Department of Biochemistry, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland.
| | - Amedeo Caflisch
- University of Zurich, Department of Biochemistry, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland.
| |
Collapse
|