1
|
Tamagnone S, Laio A, Gabrié M. Coarse-Grained Molecular Dynamics with Normalizing Flows. J Chem Theory Comput 2024. [PMID: 39223750 DOI: 10.1021/acs.jctc.4c00700] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 09/04/2024]
Abstract
We propose a sampling algorithm relying on a collective variable (CV) of midsize dimension modeled by a normalizing flow and using nonequilibrium dynamics to propose full configurational moves from the proposition of a refreshed value of the CV made by the flow. The algorithm takes the form of a Markov chain with nonlocal updates, allowing jumps through energy barriers across metastable states. The flow is trained throughout the algorithm to reproduce the free energy landscape of the CV. The output of the algorithm is a sample of thermalized configurations and the trained network that can be used to efficiently produce more configurations. We show the functioning of the algorithm first in a test case with a mixture of Gaussians. Then, we successfully tested it on a higher-dimensional system consisting of a polymer in solution with a compact state and an extended stable state separated by a high free energy barrier.
Collapse
Affiliation(s)
- Samuel Tamagnone
- International School for Advanced Studies (SISSA), Via Bonomea 265, Trieste 34136, Italy
| | - Alessandro Laio
- International School for Advanced Studies (SISSA), Via Bonomea 265, Trieste 34136, Italy
- The Abdus Salam International Centre for Theoretical Physics (ICTP), Strada Costiera 11, Trieste 34151, Italy
| | - Marylou Gabrié
- CMAP, CNRS, Institut Polytechnique de Paris, École Polytechnique, 91120 Palaiseau, France
| |
Collapse
|
2
|
Galliano L, Rende R, Coslovich D. Policy-guided Monte Carlo on general state spaces: Application to glass-forming mixtures. J Chem Phys 2024; 161:064503. [PMID: 39132794 DOI: 10.1063/5.0221221] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 05/30/2024] [Accepted: 07/23/2024] [Indexed: 08/13/2024] Open
Abstract
Policy-guided Monte Carlo is an adaptive method to simulate classical interacting systems. It adjusts the proposal distribution of the Metropolis-Hastings algorithm to maximize the sampling efficiency, using a formalism inspired by reinforcement learning. In this work, we first extend the policy-guided method to deal with a general state space, comprising, for instance, both discrete and continuous degrees of freedom, and then apply it to a few paradigmatic models of glass-forming mixtures. We assess the efficiency of a set of physically inspired moves whose proposal distributions are optimized through on-policy learning. Compared to conventional Monte Carlo methods, the optimized proposals are two orders of magnitude faster for an additive soft sphere mixture but yield a much more limited speed-up for the well-studied Kob-Andersen model. We discuss the current limitations of the method and suggest possible ways to improve it.
Collapse
Affiliation(s)
- Leonardo Galliano
- Dipartimento di Fisica, Università di Trieste, Strada Costiera 11, 34151 Trieste, Italy
| | - Riccardo Rende
- International School for Advanced Studies (SISSA), Via Bonomea 265, 34136 Trieste, Italy
| | - Daniele Coslovich
- Dipartimento di Fisica, Università di Trieste, Strada Costiera 11, 34151 Trieste, Italy
| |
Collapse
|
3
|
Kim J, Rotenberg B. Donnan equilibrium in charged slit-pores from a hybrid nonequilibrium molecular dynamics/Monte Carlo method with ions and solvent exchange. J Chem Phys 2024; 161:054107. [PMID: 39087531 DOI: 10.1063/5.0220913] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 05/29/2024] [Accepted: 07/14/2024] [Indexed: 08/02/2024] Open
Abstract
Ion partitioning between different compartments (e.g., a porous material and a bulk solution reservoir), known as Donnan equilibrium, plays a fundamental role in various contexts such as energy, environment, or water treatment. The linearized Poisson-Boltzmann (PB) equation, capturing the thermal motion of the ions with mean-field electrostatic interactions, is practically useful to understand and predict ion partitioning, despite its limited applicability to conditions of low salt concentrations and surface charge densities. Here, we investigate the Donnan equilibrium of coarse-grained dilute electrolytes confined in charged slit-pores in equilibrium with a reservoir of ions and solvent. We introduce and use an extension to confined systems of a recently developed hybrid nonequilibrium molecular dynamics/grand canonical Monte Carlo simulation method ("H4D"), which enhances the efficiency of solvent and ion-pair exchange via a fourth spatial dimension. We show that the validity range of linearized PB theory to predict the Donnan equilibrium of dilute electrolytes can be extended to highly charged pores by simply considering renormalized surface charge densities. We compare with simulations of implicit solvent models of electrolytes and show that in the low salt concentrations and thin electric double layer limit considered here, an explicit solvent has a limited effect on the Donnan equilibrium and that the main limitations of the analytical predictions are not due to the breakdown of the mean-field description but rather to the charge renormalization approximation, because it only focuses on the behavior far from the surfaces.
Collapse
Affiliation(s)
- Jeongmin Kim
- Department of Energy Engineering, Korea Institute of Energy Technology (KENTECH), Naju 58330, Republic of Korea
| | - Benjamin Rotenberg
- Sorbonne Université, CNRS, Physico-chimie des Électrolytes et Nanosystèmes Interfaciaux, PHENIX, F-75005 Paris, France
- Réseau sur le Stockage Electrochimique de l'Energie (RS2E), FR CNRS 3459, 80039 Amiens Cedex, France
| |
Collapse
|
4
|
Szczepaniak F, Dehez F, Roux B. Configurational Sampling of All-Atom Solvated Membranes Using Hybrid Nonequilibrium Molecular Dynamics Monte Carlo Simulations. J Phys Chem Lett 2024; 15:3796-3804. [PMID: 38557030 DOI: 10.1021/acs.jpclett.4c00305] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 04/04/2024]
Abstract
All-atom simulations are a powerful approach to study the structure and dynamics of biological membranes. However, sampling the atomic configurations of inhomogeneous membranes can be challenging due to the slow lateral diffusion of the various constituents. To address this issue, a hybrid nonequilibrium molecular dynamics Monte Carlo (neMD/MC) simulation method is proposed in which randomly chosen lipid molecules are swapped to generate configurations that are subsequently accepted or rejected according to a Metropolis criterion based on the alchemical work for the attempted swap calculated via a short trajectory. A dual-topology framework constraining the common atoms of the exchanging molecules yields a good acceptance probability using switching trajectories as short as 10 ps. The performance of the hybrid neMD/MC algorithm and its ability to sample the distribution of lipids near a transmembrane helix carrying a net charge are illustrated for a binary mixture of charged and zwitterionic lipids.
Collapse
Affiliation(s)
- Florence Szczepaniak
- CNRS, LPCT, Université de Lorraine, F-54000 Nancy, France
- Department of Chemistry, University of Chicago, Chicago, Illinois 60637,United States
| | - François Dehez
- CNRS, LPCT, Université de Lorraine, F-54000 Nancy, France
- Laboratoire International Associé Centre National de la Recherche Scientifique et University of Illinois at Urbana-Champaign, Université de Lorraine, LPCT, F-54000 Nancy, France
| | - Benoît Roux
- Department of Chemistry, University of Chicago, Chicago, Illinois 60637,United States
| |
Collapse
|
5
|
Deng J, Cui Q. Efficient Sampling of Cavity Hydration in Proteins with Nonequilibrium Grand Canonical Monte Carlo and Polarizable Force Fields. J Chem Theory Comput 2024; 20:1897-1911. [PMID: 38417108 DOI: 10.1021/acs.jctc.4c00013] [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: 03/01/2024]
Abstract
Prediction of the hydration levels of protein cavities and active sites is important to both mechanistic analysis and ligand design. Due to the unique microscopic environment of these buried water molecules, a polarizable model is expected to be crucial for an accurate treatment of protein internal hydration in simulations. Here we adapt a nonequilibrium candidate Monte Carlo approach for conducting grand canonical Monte Carlo simulations with the Drude polarizable force field. The GPU implementation enables the efficient sampling of internal cavity hydration levels in biomolecular systems. We also develop an enhanced sampling approach referred to as B-walking, which satisfies detailed balance and readily combines with grand canonical integration to efficiently calculate quantitative binding free energies of water to protein cavities. Applications of these developments are illustrated in a solvent box and the polar ligand binding site in trypsin. Our simulation results show that including electronic polarization leads to a modest but clear improvement in the description of water position and occupancy compared to the crystal structure. The B-walking approach enhances the range of water sampling in different chemical potential windows and thus improves the accuracy of water binding free energy calculations.
Collapse
Affiliation(s)
- Jiahua Deng
- Department of Chemistry, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, United States
| | - Qiang Cui
- Department of Chemistry, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, United States
- Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, United States
- Department of Biomedical Engineering, Boston University, 44 Cummington Mall, Boston, Massachusetts 02215, United States
| |
Collapse
|
6
|
Liu K, Rotskoff GM, Vanden-Eijnden E, Hocky GM. Computing equilibrium free energies through a nonequilibrium quench. J Chem Phys 2024; 160:034109. [PMID: 38240301 PMCID: PMC10799689 DOI: 10.1063/5.0176700] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 09/15/2023] [Accepted: 12/05/2023] [Indexed: 01/22/2024] Open
Abstract
Many methods to accelerate sampling of molecular configurations are based on the idea that temperature can be used to accelerate rare transitions. These methods typically compute equilibrium properties at a target temperature using reweighting or through Monte Carlo exchanges between replicas at higher temperatures. A recent paper [G. M. Rotskoff and E. Vanden-Eijnden, Phys. Rev. Lett. 122, 150602 (2019)] demonstrated that accurate equilibrium densities of states can also be computed through a nonequilibrium "quench" process, where sampling is performed at a higher temperature to encourage rapid mixing and then quenched to lower energy states with dissipative dynamics. Here, we provide an implementation of the quench dynamics in LAMMPS and evaluate a new formulation of nonequilibrium estimators for the computation of partition functions or free energy surfaces (FESs) of molecular systems. We show that the method is exact for a minimal model of N-independent harmonic springs and use these analytical results to develop heuristics for the amount of quenching required to obtain accurate sampling. We then test the quench approach on alanine dipeptide, where we show that it gives an FES that is accurate near the most stable configurations using the quench approach but disagrees with a reference umbrella sampling calculation in high FE regions. We then show that combining quenching with umbrella sampling allows the efficient calculation of the free energy in all regions. Moreover, by using this combined scheme, we obtain the FES across a range of temperatures at no additional cost, making it much more efficient than standard umbrella sampling if this information is required. Finally, we discuss how this approach can be extended to solute tempering and demonstrate that it is highly accurate for the case of solvated alanine dipeptide without any additional modifications.
Collapse
Affiliation(s)
| | - Grant M. Rotskoff
- Department of Chemistry, Stanford University, Stanford, California 94305, USA
| | - Eric Vanden-Eijnden
- Courant Institute of Mathematical Sciences, New York University, New York, New York 10012, USA
| | | |
Collapse
|
7
|
de Izarra A, Coudert FX, Fuchs AH, Boutin A. Molecular Simulation of the Impact of Defects on Electrolyte Intrusion in Zeolites. LANGMUIR : THE ACS JOURNAL OF SURFACES AND COLLOIDS 2023; 39:19056-19063. [PMID: 38088342 DOI: 10.1021/acs.langmuir.3c03306] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Subscribe] [Scholar Register] [Indexed: 12/27/2023]
Abstract
We have investigated through molecular simulation the intrusion of electrolytes in two representative pure-silica zeolites, silicalite-1 and chabazite, in which point defects were introduced in varying amounts. We distinguish between two types of defects, considering either "weak" or "strong" silanol nest defects, resulting in different hydration behaviors. In the presence of weak defects, the hydration process occurs through a homogeneous nucleation process, while with strong defects, we observe an initial adsorption followed by a filling of the nanoporous volume at a higher pressure. However, we show that electrolytes do not penetrate the zeolites, and these defects appear to have only marginal influence on the thermodynamics of electrolyte intrusion. While replacing pure water by the electrolyte solution shifts the intrusion pressure toward higher values because of the drop of water saturation vapor pressure, an increase in hydrophilicity of the framework due to point defects has the opposite effect, showing that controlling the amount of defects in zeolites is crucial for storage energy applications.
Collapse
Affiliation(s)
- Ambroise de Izarra
- PASTEUR, Département de Chimie, École Normale Supérieure, PSL University, Sorbonne Université, CNRS, 75005 Paris, France
- Chimie ParisTech, PSL University, CNRS, Institut de Recherche de Chimie Paris, 75005 Paris, France
| | - François-Xavier Coudert
- Chimie ParisTech, PSL University, CNRS, Institut de Recherche de Chimie Paris, 75005 Paris, France
| | - Alain H Fuchs
- Chimie ParisTech, PSL University, CNRS, Institut de Recherche de Chimie Paris, 75005 Paris, France
| | - Anne Boutin
- PASTEUR, Département de Chimie, École Normale Supérieure, PSL University, Sorbonne Université, CNRS, 75005 Paris, France
| |
Collapse
|
8
|
Kim J, Belloni L, Rotenberg B. Grand-canonical molecular dynamics simulations powered by a hybrid 4D nonequilibrium MD/MC method: Implementation in LAMMPS and applications to electrolyte solutions. J Chem Phys 2023; 159:144802. [PMID: 37819001 DOI: 10.1063/5.0168878] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 07/21/2023] [Accepted: 09/19/2023] [Indexed: 10/13/2023] Open
Abstract
Molecular simulations in an open environment, involving ion exchange, are necessary to study various systems, from biosystems to confined electrolytes. However, grand-canonical simulations are often computationally demanding in condensed phases. A promising method [L. Belloni, J. Chem. Phys. 151, 021101 (2019)], one of the hybrid nonequilibrium molecular dynamics/Monte Carlo algorithms, was recently developed, which enables efficient computation of fluctuating number or charge density in dense fluids or ionic solutions. This method facilitates the exchange through an auxiliary dimension, orthogonal to all physical dimensions, by reducing initial steric and electrostatic clashes in three-dimensional systems. Here, we report the implementation of the method in LAMMPS with a Python interface, allowing facile access to grand-canonical molecular dynamics simulations with massively parallelized computation. We validate our implementation with two electrolytes, including a model Lennard-Jones electrolyte similar to a restricted primitive model and aqueous solutions. We find that electrostatic interactions play a crucial role in the overall efficiency due to their long-range nature, particularly for water or ion-pair exchange in aqueous solutions. With properly screened electrostatic interactions and bias-based methods, our approach enhances the efficiency of salt-pair exchange in Lennard-Jones electrolytes by approximately four orders of magnitude, compared to conventional grand-canonical Monte Carlo. Furthermore, the acceptance rate of NaCl-pair exchange in aqueous solutions at moderate concentrations reaches about 3% at the maximum efficiency.
Collapse
Affiliation(s)
- Jeongmin Kim
- Sorbonne Université, CNRS, Physico-Chimie des Électrolytes et Nanosystèmes Interfaciaux, PHENIX, F-75005 Paris, France
| | - Luc Belloni
- LIONS, NIMBE, CEA, CNRS, Université Paris-Saclay, 91191 Gif-sur-Yvette, France
| | - Benjamin Rotenberg
- Sorbonne Université, CNRS, Physico-Chimie des Électrolytes et Nanosystèmes Interfaciaux, PHENIX, F-75005 Paris, France
- Réseau sur le Stockage Électrochimique de Énergie (RS2E), FR CNRS 3459, 80039 Amiens Cedex, France
| |
Collapse
|
9
|
Willow SY, Kang L, Minh DDL. Learned mappings for targeted free energy perturbation between peptide conformations. J Chem Phys 2023; 159:124104. [PMID: 38127367 PMCID: PMC10517865 DOI: 10.1063/5.0164662] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 06/23/2023] [Accepted: 09/04/2023] [Indexed: 12/23/2023] Open
Abstract
Targeted free energy perturbation uses an invertible mapping to promote configuration space overlap and the convergence of free energy estimates. However, developing suitable mappings can be challenging. Wirnsberger et al. [J. Chem. Phys. 153, 144112 (2020)] demonstrated the use of machine learning to train deep neural networks that map between Boltzmann distributions for different thermodynamic states. Here, we adapt their approach to the free energy differences of a flexible bonded molecule, deca-alanine, with harmonic biases and different spring centers. When the neural network is trained until "early stopping"-when the loss value of the test set increases-we calculate accurate free energy differences between thermodynamic states with spring centers separated by 1 Å and sometimes 2 Å. For more distant thermodynamic states, the mapping does not produce structures representative of the target state, and the method does not reproduce reference calculations.
Collapse
Affiliation(s)
- Soohaeng Yoo Willow
- Department of Chemistry, Illinois Institute of Technology, Chicago, Illinois 60616, USA
| | - Lulu Kang
- Department of Applied Mathematics, Illinois Institute of Technology, Chicago, Illinois 60616, USA
| | - David D. L. Minh
- Department of Chemistry, Department of Biology, and Center for Interdisciplinary Scientific Computation, Illinois Institute of Technology, Chicago, Illinois 60616, USA
| |
Collapse
|
10
|
Zhao M, Kognole AA, Jo S, Tao A, Hazel A, MacKerell AD. GPU-specific algorithms for improved solute sampling in grand canonical Monte Carlo simulations. J Comput Chem 2023; 44:1719-1732. [PMID: 37093676 PMCID: PMC10330275 DOI: 10.1002/jcc.27121] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Key Words] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 11/14/2022] [Revised: 03/22/2023] [Accepted: 04/12/2023] [Indexed: 04/25/2023]
Abstract
The Grand Canonical Monte Carlo (GCMC) ensemble defined by the excess chemical potential, μex , volume, and temperature, in the context of molecular simulations allows for variations in the number of particles in the system. In practice, GCMC simulations have been widely applied for the sampling of rare gasses and water, but limited in the context of larger molecules. To overcome this limitation, the oscillating μex GCMC method was introduced and shown to be of utility for sampling small solutes, such as formamide, propane, and benzene, as well as for ionic species such as monocations, acetate, and methylammonium. However, the acceptance of GCMC insertions is low, and the method is computationally demanding. In the present study, we improved the sampling efficiency of the GCMC method using known cavity-bias and configurational-bias algorithms in the context of GPU architecture. Specifically, for GCMC simulations of aqueous solution systems, the configurational-bias algorithm was extended by applying system partitioning in conjunction with a random interval extraction algorithm, thereby improving the efficiency in a highly parallel computing environment. The method is parallelized on the GPU using CUDA and OpenCL, allowing for the code to run on both Nvidia and AMD GPUs, respectively. Notably, the method is particularly well suited for GPU computing as the large number of threads allows for simultaneous sampling of a large number of configurations during insertion attempts without additional computational overhead. In addition, the partitioning scheme allows for simultaneous insertion attempts for large systems, offering considerable efficiency. Calculations on the BK Channel, a transporter, including a lipid bilayer with over 760,000 atoms, show a speed up of ~53-fold through the use of system partitioning. The improved algorithm is then combined with an enhanced μex oscillation protocol and shown to be of utility in the context of the site-identification by ligand competitive saturation (SILCS) co-solvent sampling approach as illustrated through application to the protein CDK2.
Collapse
Affiliation(s)
- Mingtian Zhao
- Computer Aided Drug Design Center, Department of Pharmaceutical Sciences, University of Maryland School of Pharmacy, Baltimore, Maryland, USA
| | | | | | | | - Anthony Hazel
- Computer Aided Drug Design Center, Department of Pharmaceutical Sciences, University of Maryland School of Pharmacy, Baltimore, Maryland, USA
| | - Alexander D MacKerell
- Computer Aided Drug Design Center, Department of Pharmaceutical Sciences, University of Maryland School of Pharmacy, Baltimore, Maryland, USA
| |
Collapse
|
11
|
Gracia Carmona O, Gillhofer M, Tomasiak L, De Ruiter A, Oostenbrink C. Accelerated Enveloping Distribution Sampling to Probe the Presence of Water Molecules. J Chem Theory Comput 2023. [PMID: 37167545 DOI: 10.1021/acs.jctc.3c00109] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 05/13/2023]
Abstract
Determining the presence of water molecules at protein-ligand interfaces is still a challenging task in free-energy calculations. The inappropriate placement of water molecules results in the stabilization of wrong conformational orientations of the ligand. With classical alchemical perturbation methods, such as thermodynamic integration (TI), it is essential to know the amount of water molecules in the active site of the respective ligands. However, the resolution of the crystal structure and the correct assignment of the electron density do not always lead to a clear placement of water molecules. In this work, we apply the one-step perturbation method named accelerated enveloping distribution sampling (AEDS) to determine the presence of water molecules in the active site by probing them in a fast and straightforward way. Based on these results, we combined the AEDS method with standard TI to calculate accurate binding free energies in the presence of buried water molecules. The main idea is to perturb the water molecules with AEDS such that they are allowed to alternate between regular water molecules and non-interacting dummy particles while treating the ligand with TI over an alchemical pathway. We demonstrate the use of AEDS to probe the presence of water molecules for six different test systems. For one of these, previous calculations showed difficulties to reproduce the experimental binding free energies, and here, we use the combined TI-AEDS approach to tackle these issues.
Collapse
Affiliation(s)
- Oriol Gracia Carmona
- Institute for Molecular Modeling and Simulation, Department of Material Sciences and Process Engineering, University of Natural Resources and Life Sciences, Vienna, Muthgasse 18, 1190 Vienna, Austria
| | - Michael Gillhofer
- Institute for Molecular Modeling and Simulation, Department of Material Sciences and Process Engineering, University of Natural Resources and Life Sciences, Vienna, Muthgasse 18, 1190 Vienna, Austria
| | - Lisa Tomasiak
- Institute for Molecular Modeling and Simulation, Department of Material Sciences and Process Engineering, University of Natural Resources and Life Sciences, Vienna, Muthgasse 18, 1190 Vienna, Austria
| | - Anita De Ruiter
- Institute for Molecular Modeling and Simulation, Department of Material Sciences and Process Engineering, University of Natural Resources and Life Sciences, Vienna, Muthgasse 18, 1190 Vienna, Austria
| | - Chris Oostenbrink
- Institute for Molecular Modeling and Simulation, Department of Material Sciences and Process Engineering, University of Natural Resources and Life Sciences, Vienna, Muthgasse 18, 1190 Vienna, Austria
- Christian Doppler Laboratory for Molecular Informatics in the Biosciences, University of Natural Resources and Life Sciences, Vienna, Muthgasse 18, 1190 Vienna, Austria
| |
Collapse
|
12
|
Melling O, Samways ML, Ge Y, Mobley DL, Essex JW. Enhanced Grand Canonical Sampling of Occluded Water Sites Using Nonequilibrium Candidate Monte Carlo. J Chem Theory Comput 2023; 19:1050-1062. [PMID: 36692215 PMCID: PMC9933432 DOI: 10.1021/acs.jctc.2c00823] [Citation(s) in RCA: 3] [Impact Index Per Article: 3.0] [Reference Citation Analysis] [Abstract] [Grants] [Track Full Text] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 08/10/2022] [Indexed: 01/25/2023]
Abstract
Water molecules play a key role in many biomolecular systems, particularly when bound at protein-ligand interfaces. However, molecular simulation studies on such systems are hampered by the relatively long time scales over which water exchange between a protein and solvent takes place. Grand canonical Monte Carlo (GCMC) is a simulation technique that avoids this issue by attempting the insertion and deletion of water molecules within a given structure. The approach is constrained by low acceptance probabilities for insertions in congested systems, however. To address this issue, here, we combine GCMC with nonequilibium candidate Monte Carlo (NCMC) to yield a method that we refer to as grand canonical nonequilibrium candidate Monte Carlo (GCNCMC), in which the water insertions and deletions are carried out in a gradual, nonequilibrium fashion. We validate this new approach by comparing GCNCMC and GCMC simulations of bulk water and three protein binding sites. We find that not only is the efficiency of the water sampling improved by GCNCMC but that it also results in increased sampling of ligand conformations in a protein binding site, revealing new water-mediated ligand-binding geometries that are not observed using alternative enhanced sampling techniques.
Collapse
Affiliation(s)
- Oliver
J. Melling
- School
of Chemistry, University of Southampton, SouthamptonSO17 1BJ, U.K.
| | - Marley L. Samways
- School
of Chemistry, University of Southampton, SouthamptonSO17 1BJ, U.K.
| | - Yunhui Ge
- Department
of Pharmaceutical Sciences, University of
California, Irvine, California92697, United States
| | - David L. Mobley
- Department
of Pharmaceutical Sciences, University of
California, Irvine, California92697, United States
- Department
of Chemistry, University of California, Irvine, California92697, United States
| | - Jonathan W. Essex
- School
of Chemistry, University of Southampton, SouthamptonSO17 1BJ, U.K.
| |
Collapse
|
13
|
Izarra AD, Coudert FX, Fuchs AH, Boutin A. Alchemical Osmostat for Monte Carlo Simulation: Sampling Aqueous Electrolyte Solution in Open Systems. J Phys Chem B 2023; 127:766-776. [PMID: 36634303 DOI: 10.1021/acs.jpcb.2c07902] [Citation(s) in RCA: 2] [Impact Index Per Article: 2.0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/14/2023]
Abstract
Molecular simulations involving electrolytes are usually performed at a fixed amount of salt ions in the simulation box, reproducing macroscopic concentration. Although this statement is valid in the bulk, the concentration of an electrolyte confined in nanoporous materials such as MOFs or zeolites is greatly affected and remains a priori unknown. The nanoporous material in equilibrium with the bulk electrolyte exchange water and ions at a given chemical potential Δμ in the semi-grand-canonical ensemble, that must be calibrated in order to determine the concentration in the nanoporous material. In this work, we propose an algorithm based on nonequilibrium candidate Monte Carlo (NCMC) moves to ultimately perform MC simulations in contact with a saline reservoir. First, we adapt the Widom insertion technique to calibrate the chemical potential by alchemically transmuting water molecules into ions by using NCMC moves. The chemical potential defines a Monte Carlo osmostat in the semi-grand-constant volume and temperature ensemble (Δμ, N, V, T) to be added in a Monte Carlo simulation where the number of ions fluctuates. In order to validate the method, we adapted the NCMC move to determine the free energy of water solvation and subsequently explore thermodynamics of electrolyte solvation at infinite dilution in water. Finally, we implemented the osmostat in MC simulations initialized with bulk water that are driven toward electrolyte solutions of similar concentration as the saline reservoir. Our results demonstrate that alchemical osmostat for MC simulation is a promising tool for use to sample electrolyte insertion in nanoporous materials.
Collapse
Affiliation(s)
- Ambroise de Izarra
- PASTEUR, Département de Chimie, École Normale Supérieure, PSL University, Sorbonne Université, CNRS, 75005, Paris, France.,Chimie ParisTech, PSL University, CNRS, Institut de Recherche de Chimie Paris, Paris75005, France
| | - François-Xavier Coudert
- Chimie ParisTech, PSL University, CNRS, Institut de Recherche de Chimie Paris, Paris75005, France
| | - Alain H Fuchs
- Chimie ParisTech, PSL University, CNRS, Institut de Recherche de Chimie Paris, Paris75005, France
| | - Anne Boutin
- PASTEUR, Département de Chimie, École Normale Supérieure, PSL University, Sorbonne Université, CNRS, 75005, Paris, France
| |
Collapse
|
14
|
Invernizzi M, Krämer A, Clementi C, Noé F. Skipping the Replica Exchange Ladder with Normalizing Flows. J Phys Chem Lett 2022; 13:11643-11649. [PMID: 36484770 DOI: 10.1021/acs.jpclett.2c03327] [Citation(s) in RCA: 4] [Impact Index Per Article: 2.0] [Reference Citation Analysis] [Abstract] [MESH Headings] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 06/17/2023]
Abstract
We combine replica exchange (parallel tempering) with normalizing flows, a class of deep generative models. These two sampling strategies complement each other, resulting in an efficient method for sampling molecular systems characterized by rare events, which we call learned replica exchange (LREX). In LREX, a normalizing flow is trained to map the configurations of the fastest-mixing replica into configurations belonging to the target distribution, allowing direct exchanges between the two without the need to simulate intermediate replicas. This can significantly reduce the computational cost compared to standard replica exchange. The proposed method also offers several advantages with respect to Boltzmann generators that directly use normalizing flows to sample the target distribution. We apply LREX to some prototypical molecular dynamics systems, highlighting the improvements over previous methods.
Collapse
Affiliation(s)
- Michele Invernizzi
- Department of Mathematics and Computer Science, Freie Universität Berlin, 14195Berlin, Germany
| | - Andreas Krämer
- Department of Mathematics and Computer Science, Freie Universität Berlin, 14195Berlin, Germany
| | - Cecilia Clementi
- Department of Physics, Freie Universität Berlin, 14195Berlin, Germany
- Department of Chemistry, Rice University, 77005Houston, United States
- Center for Theoretical Biological Physics, Rice University, 77005Houston, United States
| | - Frank Noé
- Department of Mathematics and Computer Science, Freie Universität Berlin, 14195Berlin, Germany
- Department of Physics, Freie Universität Berlin, 14195Berlin, Germany
- Department of Chemistry, Rice University, 77005Houston, United States
- Microsoft Research AI4Science, 10178Berlin, Germany
| |
Collapse
|
15
|
Abbas G, Cardenas AE, Elber R. The Structures of Heterogeneous Membranes and Their Interactions with an Anticancer Peptide: A Molecular Dynamics Study. Life (Basel) 2022; 12:1473. [PMID: 36294908 PMCID: PMC9604715 DOI: 10.3390/life12101473] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Key Words] [Grants] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 08/31/2022] [Revised: 09/16/2022] [Accepted: 09/20/2022] [Indexed: 12/02/2022] Open
Abstract
We conduct molecular dynamics simulations of model heterogeneous membranes and their interactions with a 24-amino acid peptide-NAF-144-67. NAF-144-67 is an anticancer peptide that selectively permeates and kills malignant cells; it does not permeate normal cells. We examine three membranes with different binary mixtures of lipids, DOPC-DOPA, DOPC-DOPS, and DOPC-DOPE, with a single peptide embedded in each as models for the diversity of biological membranes. We illustrate that the peptide organization in the membrane depends on the types of nearby phospholipids and is influenced by the charge and size of the head groups. The present study sheds light on early events of permeation and the mechanisms by which an amphiphilic peptide crosses from an aqueous solution to a hydrophobic membrane. Understanding the translocation mechanism is likely to help the design of new permeants.
Collapse
Affiliation(s)
- Ghulam Abbas
- National Center for Bioinformatics, Quaid-i-Azam University, Islamabad 45320, Pakistan or
| | - Alfredo E. Cardenas
- Oden Institute for Computational and Engineering Sciences, University of Texas at Austin, Austin, TX 78712, USA
| | - Ron Elber
- Oden Institute for Computational and Engineering Sciences, University of Texas at Austin, Austin, TX 78712, USA
- Department of Chemistry, University of Texas at Austin, Austin, TX 78712, USA
| |
Collapse
|
16
|
Das P, Parmar ADS, Sastry S. Annealing glasses by cyclic shear deformation. J Chem Phys 2022; 157:044501. [DOI: 10.1063/5.0100523] [Citation(s) in RCA: 2] [Impact Index Per Article: 1.0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/24/2023] Open
Abstract
A major challenge in simulating glassy systems is the ability to generate configurations that may be found in equilibrium at sufficiently low temperatures, in order to probe static and dynamic behavior close to the glass transition. A variety of approaches have recently explored ways of surmounting this obstacle. Here, we explore the possibility of employing mechanical agitation, in the form of cyclic shear deformation, to generate low energy configurations in a model glass former. We perform shear deformation simulations over a range of temperatures, shear rates, and strain amplitudes. We find that shear deformation induces faster relaxation toward low energy configurations, or overaging, in simulations at sufficiently low temperatures, consistently with previous results for athermal shear. However, for temperatures at which simulations can be run until a steady state is reached with or without shear deformation, we find that the inclusion of shear deformation does not result in any speed up of the relaxation toward low energy configurations. Although we find the configurations from shear simulations to have properties indistinguishable from an equilibrium ensemble, the cyclic shear procedure does not guarantee that we generate an equilibrium ensemble at a desired temperature. In order to ensure equilibrium sampling, we develop a hybrid Monte Carlo algorithm that employs cyclic shear as a trial generation step and has acceptance probabilities that depend not only on the change in internal energy but also on the heat dissipated (equivalently, work done). We show that such an algorithm, indeed, generates an equilibrium ensemble.
Collapse
Affiliation(s)
- Pallabi Das
- Theoretical Sciences Unit and School of Advanced Materials, Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur Campus, Bengaluru 560064, India
| | - Anshul D. S. Parmar
- Theoretical Sciences Unit and School of Advanced Materials, Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur Campus, Bengaluru 560064, India
| | - Srikanth Sastry
- Theoretical Sciences Unit and School of Advanced Materials, Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur Campus, Bengaluru 560064, India
| |
Collapse
|
17
|
Patel LA, Chau P, Debesai S, Darwin L, Neale C. Drug Discovery by Automated Adaptation of Chemical Structure and Identity. J Chem Theory Comput 2022; 18:5006-5024. [PMID: 35834740 DOI: 10.1021/acs.jctc.1c01271] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/29/2022]
Abstract
Computer-aided drug design offers the potential to dramatically reduce the cost and effort required for drug discovery. While screening-based methods are valuable in the early stages of hit identification, they are frequently succeeded by iterative, hypothesis-driven computations that require recurrent investment of human time and intuition. To increase automation, we introduce a computational method for lead refinement that combines concerted dynamics of the ligand/protein complex via molecular dynamics simulations with integrated Monte Carlo-based changes in the chemical formula of the ligand. This approach, which we refer to as ligand-exchange Monte Carlo molecular dynamics, accounts for solvent- and entropy-based contributions to competitive binding free energies by coupling the energetics of bound and unbound states during the ligand-exchange attempt. Quantitative comparison of relative binding free energies to reference values from free energy perturbation, conducted in vacuum, indicates that ligand-exchange Monte Carlo molecular dynamics simulations sample relevant conformational ensembles and are capable of identifying strongly binding compounds. Additional simulations demonstrate the use of an implicit solvent model. We speculate that the use of chemical graphs in which exchanges are only permitted between ligands with sufficient similarity may enable an automated search to capture some of the benefits provided by human intuition during hypothesis-guided lead refinement.
Collapse
|
18
|
Dai C, Heng J, Jacob PE, Whiteley N. An Invitation to Sequential Monte Carlo Samplers. J Am Stat Assoc 2022. [DOI: 10.1080/01621459.2022.2087659] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 10/18/2022]
Affiliation(s)
| | | | | | - Nick Whiteley
- School of Mathematics, University of Bristol, Bristol, UK
| |
Collapse
|
19
|
Monroe JI, Shen VK. Learning Efficient, Collective Monte Carlo Moves with Variational Autoencoders. J Chem Theory Comput 2022; 18:3622-3636. [PMID: 35613327 PMCID: PMC11210279 DOI: 10.1021/acs.jctc.2c00110] [Citation(s) in RCA: 2] [Impact Index Per Article: 1.0] [Reference Citation Analysis] [Abstract] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/28/2022]
Abstract
Discovering meaningful collective variables for enhancing sampling, via applied biasing potentials or tailored MC move sets, remains a major challenge within molecular simulation. While recent studies identifying collective variables with variational autoencoders (VAEs) have focused on the encoding and latent space discovered by a VAE, the impact of the decoding and its ability to act as a generative model remains unexplored. We demonstrate how VAEs may be used to learn (on-the-fly and with minimal human intervention) highly efficient, collective Monte Carlo moves that accelerate sampling along the learned collective variable. In contrast to many machine learning-based efforts to bias sampling and generate novel configurations, our methods result in exact sampling in the ensemble of interest and do not require reweighting. In fact, we show that the acceptance rates of our moves approach unity for a perfect VAE model. While this is never observed in practice, VAE-based Monte Carlo moves still enhance sampling of new configurations. We demonstrate, however, that the form of the encoding and decoding distributions, in particular the extent to which the decoder reflects the underlying physics, greatly impacts the performance of the trained VAE.
Collapse
Affiliation(s)
- Jacob I Monroe
- Chemical Sciences Division, National Institute of Standards and Technology, Gaithersburg, Maryland 20899-8320, United States
| | - Vincent K Shen
- Chemical Sciences Division, National Institute of Standards and Technology, Gaithersburg, Maryland 20899-8320, United States
| |
Collapse
|
20
|
Suruzhon M, Bodnarchuk MS, Ciancetta A, Wall ID, Essex JW. Enhancing Ligand and Protein Sampling Using Sequential Monte Carlo. J Chem Theory Comput 2022; 18:3894-3910. [PMID: 35588256 PMCID: PMC9202307 DOI: 10.1021/acs.jctc.1c01198] [Citation(s) in RCA: 0] [Impact Index Per Article: 0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/28/2022]
Abstract
The sampling problem is one of the most widely studied topics in computational chemistry. While various methods exist for sampling along a set of reaction coordinates, many require system-dependent hyperparameters to achieve maximum efficiency. In this work, we present an alchemical variation of adaptive sequential Monte Carlo (SMC), an irreversible importance resampling method that is part of a well-studied class of methods that have been used in various applications but have been underexplored in computational biophysics. Afterward, we apply alchemical SMC on a variety of test cases, including torsional rotations of solvated ligands (butene and a terphenyl derivative), translational and rotational movements of protein-bound ligands, and protein side chain rotation coupled to the ligand degrees of freedom (T4-lysozyme, protein tyrosine phosphatase 1B, and transforming growth factor β). We find that alchemical SMC is an efficient way to explore targeted degrees of freedom and can be applied to a variety of systems using the same hyperparameters to achieve a similar performance. Alchemical SMC is a promising tool for preparatory exploration of systems where long-timescale sampling of the entire system can be traded off against short-timescale sampling of a particular set of degrees of freedom over a population of conformers.
Collapse
Affiliation(s)
- Miroslav Suruzhon
- School of Chemistry, University of Southampton, Highfield, Southampton SO17 1BJ, U.K
| | | | - Antonella Ciancetta
- Sygnature Discovery, Bio City, Pennyfoot Street, Nottingham NG1 1GR, U.K.,Department of Chemical, Pharmaceutical and Agricultural Sciences─DOCPAS, University of Ferrara, Via Fossato di Mortara 17/19, 44121 Ferrara, Italy
| | - Ian D Wall
- GSK Medicines Research Centre, Gunnels Wood Road, Stevenage SG1 2NY, U.K
| | - Jonathan W Essex
- School of Chemistry, University of Southampton, Highfield, Southampton SO17 1BJ, U.K
| |
Collapse
|
21
|
Wang S, Venkatesh A, Ramkrishna D, Narsimhan V. Brownian bridges for stochastic chemical processes-An approximation method based on the asymptotic behavior of the backward Fokker-Planck equation. J Chem Phys 2022; 156:184108. [PMID: 35568530 DOI: 10.1063/5.0080540] [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
A Brownian bridge is a continuous random walk conditioned to end in a given region by adding an effective drift to guide paths toward the desired region of phase space. This idea has many applications in chemical science where one wants to control the endpoint of a stochastic process-e.g., polymer physics, chemical reaction pathways, heat/mass transfer, and Brownian dynamics simulations. Despite its broad applicability, the biggest limitation of the Brownian bridge technique is that it is often difficult to determine the effective drift as it comes from a solution of a Backward Fokker-Planck (BFP) equation that is infeasible to compute for complex or high-dimensional systems. This paper introduces a fast approximation method to generate a Brownian bridge process without solving the BFP equation explicitly. Specifically, this paper uses the asymptotic properties of the BFP equation to generate an approximate drift and determine ways to correct (i.e., re-weight) any errors incurred from this approximation. Because such a procedure avoids the solution of the BFP equation, we show that it drastically accelerates the generation of conditioned random walks. We also show that this approach offers reasonable improvement compared to other sampling approaches using simple bias potentials.
Collapse
Affiliation(s)
- Shiyan Wang
- Davidson School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907, USA
| | - Anirudh Venkatesh
- Davidson School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907, USA
| | - Doraiswami Ramkrishna
- Davidson School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907, USA
| | - Vivek Narsimhan
- Davidson School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907, USA
| |
Collapse
|
22
|
Ge Y, Wych DC, Samways ML, Wall ME, Essex JW, Mobley DL. Enhancing Sampling of Water Rehydration on Ligand Binding: A Comparison of Techniques. J Chem Theory Comput 2022; 18:1359-1381. [PMID: 35148093 PMCID: PMC9241631 DOI: 10.1021/acs.jctc.1c00590] [Citation(s) in RCA: 20] [Impact Index Per Article: 10.0] [Reference Citation Analysis] [Abstract] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/18/2023]
Abstract
Water often plays a key role in protein structure, molecular recognition, and mediating protein-ligand interactions. Thus, free energy calculations must adequately sample water motions, which often proves challenging in typical MD simulation time scales. Thus, the accuracy of methods relying on MD simulations ends up limited by slow water sampling. Particularly, as a ligand is removed or modified, bulk water may not have time to fill or rearrange in the binding site. In this work, we focus on several molecular dynamics (MD) simulation-based methods attempting to help rehydrate buried water sites: BLUES, using nonequilibrium candidate Monte Carlo (NCMC); grand, using grand canonical Monte Carlo (GCMC); and normal MD. We assess the accuracy and efficiency of these methods in rehydrating target water sites. We selected a range of systems with varying numbers of waters in the binding site, as well as those where water occupancy is coupled to the identity or binding mode of the ligand. We analyzed the rehydration of buried water sites in binding pockets using both clustering of trajectories and direct analysis of electron density maps. Our results suggest both BLUES and grand enhance water sampling relative to normal MD and grand is more robust than BLUES, but also that water sampling remains a major challenge for all of the methods tested. The lessons we learned for these methods and systems are discussed.
Collapse
Affiliation(s)
- Yunhui Ge
- Department of Pharmaceutical Sciences, University of California, Irvine, California 92697, United States
| | - David C Wych
- Department of Pharmaceutical Sciences, University of California, Irvine, California 92697, United States
- Computer, Computational, and Statistical Sciences Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, United States
| | - Marley L Samways
- School of Chemistry, University of Southampton, Southampton SO17 1BJ, United Kingdom
| | - Michael E Wall
- Computer, Computational, and Statistical Sciences Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, United States
| | - Jonathan W Essex
- School of Chemistry, University of Southampton, Southampton SO17 1BJ, United Kingdom
| | - David L Mobley
- Department of Pharmaceutical Sciences, University of California, Irvine, California 92697, United States
- Department of Chemistry, University of California, Irvine, California 92697, United States
| |
Collapse
|
23
|
Madin OC, Boothroyd S, Messerly RA, Fass J, Chodera JD, Shirts MR. Bayesian-Inference-Driven Model Parametrization and Model Selection for 2CLJQ Fluid Models. J Chem Inf Model 2022; 62:874-889. [PMID: 35129974 PMCID: PMC9217127 DOI: 10.1021/acs.jcim.1c00829] [Citation(s) in RCA: 3] [Impact Index Per Article: 1.5] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/28/2022]
Abstract
A high level of physical detail in a molecular model improves its ability to perform high accuracy simulations but can also significantly affect its complexity and computational cost. In some situations, it is worthwhile to add complexity to a model to capture properties of interest; in others, additional complexity is unnecessary and can make simulations computationally infeasible. In this work, we demonstrate the use of Bayesian inference for molecular model selection, using Monte Carlo sampling techniques accelerated with surrogate modeling to evaluate the Bayes factor evidence for different levels of complexity in the two-centered Lennard-Jones + quadrupole (2CLJQ) fluid model. Examining three nested levels of model complexity, we demonstrate that the use of variable quadrupole and bond length parameters in this model framework is justified only for some chemistries. Through this process, we also get detailed information about the distributions and correlation of parameter values, enabling improved parametrization and parameter analysis. We also show how the choice of parameter priors, which encode previous model knowledge, can have substantial effects on the selection of models, penalizing careless introduction of additional complexity. We detail the computational techniques used in this analysis, providing a roadmap for future applications of molecular model selection via Bayesian inference and surrogate modeling.
Collapse
Affiliation(s)
- Owen C. Madin
- Department of Chemical & Biological Engineering, University of Colorado Boulder, Boulder, CO 80309
| | - Simon Boothroyd
- Boothroyd Scientific Consulting Ltd., 71-75 Shelton Street, London, Greater London, United Kingdom, WC2H 9JQ
| | | | - Josh Fass
- Computational & Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065
| | - John D. Chodera
- Department of Chemical & Biological Engineering, University of Colorado Boulder, Boulder, CO 80309
| | - Michael R. Shirts
- Department of Chemical & Biological Engineering, University of Colorado Boulder, Boulder, CO 80309
| |
Collapse
|
24
|
Lunansky G, Naberman J, van Borkulo CD, Chen C, Wang L, Borsboom D. Intervening on psychopathology networks: Evaluating intervention targets through simulations. Methods 2021; 204:29-37. [PMID: 34793976 DOI: 10.1016/j.ymeth.2021.11.006] [Citation(s) in RCA: 23] [Impact Index Per Article: 7.7] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 07/13/2021] [Revised: 10/30/2021] [Accepted: 11/11/2021] [Indexed: 01/16/2023] Open
Abstract
Identifying the different influences of symptoms in dynamic psychopathology models may hold promise for increasing treatment efficacy in clinical applications. Dynamic psychopathology models study the behavioral patterns of symptom networks, where symptoms mutually enforce each other. Interventions could be tailored to specific symptoms that are most effective at lowering symptom activity or that hinder the further development of psychopathology. Simulating interventions in psychopathology network models fits in a novel tradition where symptom-specific perturbations are used as in silico interventions. Here, we present the NodeIdentifyR algorithm (NIRA) to identify the projected most efficient, symptom-specific intervention target in a network model (i.e., the Ising model). We implemented NIRA in a freely available R package. The technique studies the projected effects of symptom-specific interventions by simulating data while symptom parameters (i.e., thresholds) are systematically altered. The projected effect of these interventions is defined in terms of the expected change in overall symptom activity across simulations. With this algorithm, it is possible to study (1) whether symptoms differ in their projected influence on the behavior of the symptom network and, if so, (2) which symptom has the largest projected effect in lowering or increasing overall symptom activation. As an illustration, we apply the algorithm to an empirical dataset containing Post-Traumatic Stress Disorder symptom assessments of participants who experienced the Wenchuan earthquake in 2008. The most important limitations of the method are discussed, as well as recommendations for future research, such as shifting towards modeling individual processes to validate these types of simulation-based intervention methods.
Collapse
Affiliation(s)
- Gabriela Lunansky
- Department of Psychological Methods, University of Amsterdam, Amsterdam, the Netherlands.
| | - Jasper Naberman
- Department of Psychological Methods, University of Amsterdam, Amsterdam, the Netherlands
| | - Claudia D van Borkulo
- Department of Psychological Methods, University of Amsterdam, Amsterdam, the Netherlands; Centre for Urban Mental Health, University of Amsterdam, The Netherlands
| | - Chen Chen
- Laboratory for Traumatic Stress Studies, CAS Key Laboratory of Mental Health, Institute of Psychology, Chinese Academy of Sciences, Beijing, China; Department of Psychology, University of Chinese Academy of Sciences, Beijing, China
| | - Li Wang
- Laboratory for Traumatic Stress Studies, CAS Key Laboratory of Mental Health, Institute of Psychology, Chinese Academy of Sciences, Beijing, China; Department of Psychology, University of Chinese Academy of Sciences, Beijing, China
| | - Denny Borsboom
- Department of Psychological Methods, University of Amsterdam, Amsterdam, the Netherlands
| |
Collapse
|
25
|
Cherniavskyi YK, Fathizadeh A, Elber R, Tieleman DP. Computer simulations of a heterogeneous membrane with enhanced sampling techniques. J Chem Phys 2021; 153:144110. [PMID: 33086798 PMCID: PMC7556882 DOI: 10.1063/5.0014176] [Citation(s) in RCA: 8] [Impact Index Per Article: 2.7] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/14/2022] Open
Abstract
Computational determination of the equilibrium state of heterogeneous phospholipid membranes is a significant challenge. We wish to explore the rich phase diagram of these multi-component systems. However, the diffusion and mixing times in membranes are long compared to typical time scales of computer simulations. Here, we evaluate the combination of the enhanced sampling techniques molecular dynamics with alchemical steps and Monte Carlo with molecular dynamics with a coarse-grained model of membranes (Martini) to reduce the number of steps and force evaluations that are needed to reach equilibrium. We illustrate a significant gain compared to straightforward molecular dynamics of the Martini model by factors between 3 and 10. The combination is a useful tool to enhance the study of phase separation and the formation of domains in biological membranes.
Collapse
Affiliation(s)
- Yevhen K Cherniavskyi
- Department of Biological Sciences and Centre for Molecular Simulation, University of Calgary, 2500 University Drive NW, Calgary, Alberta T2N 1N4, Canada
| | - Arman Fathizadeh
- The Oden Institute for Computational Engineering and Sciences, University of Texas at Austin, Austin, Texas 78712, USA
| | - Ron Elber
- The Oden Institute for Computational Engineering and Sciences, University of Texas at Austin, Austin, Texas 78712, USA
| | - D Peter Tieleman
- Department of Biological Sciences and Centre for Molecular Simulation, University of Calgary, 2500 University Drive NW, Calgary, Alberta T2N 1N4, Canada
| |
Collapse
|
26
|
Abstract
Markov chain Monte Carlo methods are a powerful tool for sampling equilibrium configurations in complex systems. One problem these methods often face is slow convergence over large energy barriers. In this work, we propose a novel method that increases convergence in systems composed of many metastable states. This method aims to connect metastable regions directly using generative neural networks in order to propose new configurations in the Markov chain and optimizes the acceptance probability of large jumps between modes in the configuration space. We provide a comprehensive theory as well as a training scheme for the network and demonstrate the method on example systems.
Collapse
Affiliation(s)
- Luigi Sbailò
- Freie Universität Berlin, Department of Mathematics and Computer Science, Arnimallee 6, 14195 Berlin, Germany
| | - Manuel Dibak
- Freie Universität Berlin, Department of Mathematics and Computer Science, Arnimallee 6, 14195 Berlin, Germany
| | - Frank Noé
- Freie Universität Berlin, Department of Mathematics and Computer Science, Arnimallee 6, 14195 Berlin, Germany
| |
Collapse
|
27
|
Bergazin TD, Ben-Shalom IY, Lim NM, Gill SC, Gilson MK, Mobley DL. Enhancing water sampling of buried binding sites using nonequilibrium candidate Monte Carlo. J Comput Aided Mol Des 2021; 35:167-177. [PMID: 32968887 PMCID: PMC7904576 DOI: 10.1007/s10822-020-00344-8] [Citation(s) in RCA: 12] [Impact Index Per Article: 4.0] [Reference Citation Analysis] [Abstract] [Key Words] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 06/04/2020] [Accepted: 09/16/2020] [Indexed: 11/26/2022]
Abstract
Water molecules can be found interacting with the surface and within cavities in proteins. However, water exchange between bulk and buried hydration sites can be slow compared to simulation timescales, thus leading to the inefficient sampling of the locations of water. This can pose problems for free energy calculations for computer-aided drug design. Here, we apply a hybrid method that combines nonequilibrium candidate Monte Carlo (NCMC) simulations and molecular dynamics (MD) to enhance sampling of water in specific areas of a system, such as the binding site of a protein. Our approach uses NCMC to gradually remove interactions between a selected water molecule and its environment, then translates the water to a new region, before turning the interactions back on. This approach of gradual removal of interactions, followed by a move and then reintroduction of interactions, allows the environment to relax in response to the proposed water translation, improving acceptance of moves and thereby accelerating water exchange and sampling. We validate this approach on several test systems including the ligand-bound MUP-1 and HSP90 proteins with buried crystallographic waters removed. We show that our BLUES (NCMC/MD) method enhances water sampling relative to normal MD when applied to these systems. Thus, this approach provides a strategy to improve water sampling in molecular simulations which may be useful in practical applications in drug discovery and biomolecular design.
Collapse
Affiliation(s)
| | - Ido Y Ben-Shalom
- Skaggs School of Pharmacy and Pharmaceutical Sciences, University of California, San Diego, La Jolla, CA, 92093, USA
| | - Nathan M Lim
- Department of Pharmaceutical Sciences, University of California, Irvine, CA, 92697, USA
| | - Sam C Gill
- Department of Chemistry, University of California, Irvine, Irvine, CA, 92697, USA
| | - Michael K Gilson
- Skaggs School of Pharmacy and Pharmaceutical Sciences, University of California, San Diego, La Jolla, CA, 92093, USA
| | - David L Mobley
- Department of Pharmaceutical Sciences, University of California, Irvine, CA, 92697, USA.
- Department of Chemistry, University of California, Irvine, Irvine, CA, 92697, USA.
| |
Collapse
|
28
|
Fathizadeh A, Valentine M, Baiz CR, Elber R. Phase Transition in a Heterogeneous Membrane: Atomically Detailed Picture. J Phys Chem Lett 2020; 11:5263-5267. [PMID: 32525318 PMCID: PMC7334090 DOI: 10.1021/acs.jpclett.0c01255] [Citation(s) in RCA: 5] [Impact Index Per Article: 1.3] [Reference Citation Analysis] [Abstract] [Key Words] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 05/02/2023]
Abstract
Membranes serve diverse functions in biological systems. Variations in their molecular compositions impact their physical properties and lead to rich phase behavior such as switching from the gel to fluid phase and/or separation to micro- and macrodomains with different molecular compositions. We present a combined computational and experimental study of the phase behavior of a mixed membrane of 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) and 1,2-dilauroyl-sn-glycero-3-phosphocholine (DLPC) molecules. This heterogeneous membrane changes from gel to fluid and shows separate domains as a function of temperature. Atomically detailed simulations provide microscopic information about these molecular assemblies. However, these systems are challenging for computations since approaching equilibrium necessitates exceptionally long molecular dynamics trajectories. We use the simulation method of MDAS (Molecular Dynamics with Alchemical Steps) to generate adequate statistics. Isotope-edited IR spectroscopy of the lipids was used to benchmark the simulations. Together, simulations and experiments provide insight into the structural and dynamical features of the phase diagram.
Collapse
|
29
|
Cezar HM, Canuto S, Coutinho K. DICE: A Monte Carlo Code for Molecular Simulation Including the Configurational Bias Monte Carlo Method. J Chem Inf Model 2020; 60:3472-3488. [DOI: 10.1021/acs.jcim.0c00077] [Citation(s) in RCA: 24] [Impact Index Per Article: 6.0] [Reference Citation Analysis] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/29/2022]
Affiliation(s)
- Henrique M. Cezar
- Instituto de Fisica, Universidade de Sao Paulo, 05508-090 Sao Paulo, SP, Brazil
| | - Sylvio Canuto
- Instituto de Fisica, Universidade de Sao Paulo, 05508-090 Sao Paulo, SP, Brazil
| | - Kaline Coutinho
- Instituto de Fisica, Universidade de Sao Paulo, 05508-090 Sao Paulo, SP, Brazil
| |
Collapse
|
30
|
Sasmal S, Gill SC, Lim NM, Mobley DL. Sampling Conformational Changes of Bound Ligands Using Nonequilibrium Candidate Monte Carlo and Molecular Dynamics. J Chem Theory Comput 2020; 16:1854-1865. [PMID: 32058713 PMCID: PMC7325746 DOI: 10.1021/acs.jctc.9b01066] [Citation(s) in RCA: 10] [Impact Index Per Article: 2.5] [Reference Citation Analysis] [Abstract] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/01/2023]
Abstract
Flexible ligands often have multiple binding modes or bound conformations that differ by rotation of a portion of the molecule around internal rotatable bonds. Knowledge of these binding modes is important for understanding the interactions stabilizing the ligand in the binding pocket, and other studies indicate it is important for calculating accurate binding affinities. In this work, we use a hybrid molecular dynamics (MD)/nonequilibrium candidate Monte Carlo (NCMC) method to sample the different binding modes of several flexible ligands and also to estimate the population distribution of the modes. The NCMC move proposal is divided into three parts. The flexible part of the ligand is alchemically turned off by decreasing the electrostatics and steric interactions gradually, followed by rotating the rotatable bond by a random angle and then slowly turning the ligand back on to its fully interacting state. The alchemical steps prior to and after the move proposal help the surrounding protein and water atoms in the binding pocket relax around the proposed ligand conformation and increase move acceptance rates. The protein-ligand system is propagated using classical MD in between the NCMC proposals. Using this MD/NCMC method, we were able to correctly reproduce the different binding modes of inhibitors binding to two kinase targets-c-Jun N-terminal kinase-1 and cyclin-dependent kinase 2-at a much lower computational cost compared to conventional MD and umbrella sampling. This method is available as a part of the BLUES software package.
Collapse
Affiliation(s)
- Sukanya Sasmal
- Department of Pharmaceutical Sciences, University of California, Irvine, California 92697, United States
| | - Samuel C Gill
- Department of Chemistry, University of California, Irvine, California 92697, United States
| | - Nathan M Lim
- Department of Pharmaceutical Sciences, University of California, Irvine, California 92697, United States
| | - David L Mobley
- Department of Chemistry, University of California, Irvine, California 92697, United States
- Department of Pharmaceutical Sciences, University of California, Irvine, California 92697, United States
| |
Collapse
|
31
|
Rybiński M, Möller S, Sunnåker M, Lormeau C, Stelling J. TopoFilter: a MATLAB package for mechanistic model identification in systems biology. BMC Bioinformatics 2020; 21:34. [PMID: 31996136 PMCID: PMC6990465 DOI: 10.1186/s12859-020-3343-y] [Citation(s) in RCA: 2] [Impact Index Per Article: 0.5] [Reference Citation Analysis] [Abstract] [Key Words] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 03/29/2019] [Accepted: 01/08/2020] [Indexed: 12/27/2022] Open
Abstract
Background To develop mechanistic dynamic models in systems biology, one often needs to identify all (or minimal) representations of the biological processes that are consistent with experimental data, out of a potentially large set of hypothetical mechanisms. However, a simple enumeration of all alternatives becomes quickly intractable when the number of model parameters grows. Selecting appropriate dynamic models out of a large ensemble of models, taking the uncertainty in our biological knowledge and in the experimental data into account, is therefore a key current problem in systems biology. Results The TopoFilter package addresses this problem in a heuristic and automated fashion by implementing the previously described topological filtering method for Bayesian model selection. It includes a core heuristic for searching the space of submodels of a parametrized model, coupled with a sampling-based exploration of the parameter space. Recent developments of the method allow to balance exhaustiveness and speed of the model space search, to efficiently re-sample parameters, to parallelize the search, and to use custom scoring functions. We use a theoretical example to motivate these features and then demonstrate TopoFilter’s applicability for a yeast signaling network with more than 250’000 possible model structures. Conclusions TopoFilter is a flexible software framework that makes Bayesian model selection and reduction efficient and scalable to network models of a complexity that represents contemporary problems in, for example, cell signaling. TopoFilter is open-source, available under the GPL-3.0 license at https://gitlab.com/csb.ethz/TopoFilter. It includes installation instructions, a quickstart guide, a description of all package options, and multiple examples.
Collapse
Affiliation(s)
- Mikołaj Rybiński
- Department of Biosystems Science and Engineering and SIB Swiss Institute of Bioinformatics, ETH Zurich, Mattenstr. 26, Basel, 4058, Switzerland.,ID Scientific IT Services, ETH Zurich, Zurich, 8092, Switzerland
| | - Simon Möller
- Department of Biosystems Science and Engineering and SIB Swiss Institute of Bioinformatics, ETH Zurich, Mattenstr. 26, Basel, 4058, Switzerland
| | - Mikael Sunnåker
- Department of Biosystems Science and Engineering and SIB Swiss Institute of Bioinformatics, ETH Zurich, Mattenstr. 26, Basel, 4058, Switzerland
| | - Claude Lormeau
- Department of Biosystems Science and Engineering and SIB Swiss Institute of Bioinformatics, ETH Zurich, Mattenstr. 26, Basel, 4058, Switzerland.,Life Science Zurich Ph.D. program "Systems Biology", Zurich, 8092, Switzerland
| | - Jörg Stelling
- Department of Biosystems Science and Engineering and SIB Swiss Institute of Bioinformatics, ETH Zurich, Mattenstr. 26, Basel, 4058, Switzerland.
| |
Collapse
|
32
|
Copperman J, Aristoff D, Makarov DE, Simpson G, Zuckerman DM. Transient probability currents provide upper and lower bounds on non-equilibrium steady-state currents in the Smoluchowski picture. J Chem Phys 2019; 151:174108. [PMID: 31703496 PMCID: PMC7043855 DOI: 10.1063/1.5120511] [Citation(s) in RCA: 3] [Impact Index Per Article: 0.6] [Reference Citation Analysis] [Abstract] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 07/18/2019] [Accepted: 10/14/2019] [Indexed: 01/04/2023] Open
Abstract
Probability currents are fundamental in characterizing the kinetics of nonequilibrium processes. Notably, the steady-state current Jss for a source-sink system can provide the exact mean-first-passage time (MFPT) for the transition from the source to sink. Because transient nonequilibrium behavior is quantified in some modern path sampling approaches, such as the "weighted ensemble" strategy, there is strong motivation to determine bounds on Jss-and hence on the MFPT-as the system evolves in time. Here, we show that Jss is bounded from above and below by the maximum and minimum, respectively, of the current as a function of the spatial coordinate at any time t for one-dimensional systems undergoing overdamped Langevin (i.e., Smoluchowski) dynamics and for higher-dimensional Smoluchowski systems satisfying certain assumptions when projected onto a single dimension. These bounds become tighter with time, making them of potential practical utility in a scheme for estimating Jss and the long time scale kinetics of complex systems. Conceptually, the bounds result from the fact that extrema of the transient currents relax toward the steady-state current.
Collapse
Affiliation(s)
- Jeremy Copperman
- Department of Biomedical Engineering, Oregon Health and Science University, Portland, Oregon 97239, USA
| | - David Aristoff
- Department of Mathematics, Colorado State University, Fort Collins, Colorado 80523, USA
| | - Dmitrii E Makarov
- Department of Chemistry and Oden Institute for Computational Engineering and Sciences, University of Texas, Austin, Texas 78712, USA
| | - Gideon Simpson
- Department of Mathematics, Drexel University, Philadelphia, Pennsylvania 19104, USA
| | - Daniel M Zuckerman
- Department of Biomedical Engineering, Oregon Health and Science University, Portland, Oregon 97239, USA
| |
Collapse
|
33
|
Lorenzoni A, Muccini M, Mercuri F. A Computational Predictive Approach for Controlling the Morphology of Functional Molecular Aggregates on Substrates. ADVANCED THEORY AND SIMULATIONS 2019. [DOI: 10.1002/adts.201900156] [Citation(s) in RCA: 5] [Impact Index Per Article: 1.0] [Reference Citation Analysis] [Track Full Text] [Subscribe] [Scholar Register] [Indexed: 01/08/2023]
Affiliation(s)
- Andrea Lorenzoni
- Istituto per lo Studio dei Materiali Nanostrutturati (ISMN)Consiglio Nazionale delle Ricerche (CNR) Via P. Gobetti 101 40129 Bologna Italy
| | - Michele Muccini
- Istituto per lo Studio dei Materiali Nanostrutturati (ISMN)Consiglio Nazionale delle Ricerche (CNR) Via P. Gobetti 101 40129 Bologna Italy
| | - Francesco Mercuri
- Istituto per lo Studio dei Materiali Nanostrutturati (ISMN)Consiglio Nazionale delle Ricerche (CNR) Via P. Gobetti 101 40129 Bologna Italy
| |
Collapse
|
34
|
Noé F, Olsson S, Köhler J, Wu H. Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning. Science 2019; 365:365/6457/eaaw1147. [DOI: 10.1126/science.aaw1147] [Citation(s) in RCA: 205] [Impact Index Per Article: 41.0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 11/29/2018] [Accepted: 07/19/2019] [Indexed: 02/01/2023]
Abstract
Computing equilibrium states in condensed-matter many-body systems, such as solvated proteins, is a long-standing challenge. Lacking methods for generating statistically independent equilibrium samples in “one shot,” vast computational effort is invested for simulating these systems in small steps, e.g., using molecular dynamics. Combining deep learning and statistical mechanics, we developed Boltzmann generators, which are shown to generate unbiased one-shot equilibrium samples of representative condensed-matter systems and proteins. Boltzmann generators use neural networks to learn a coordinate transformation of the complex configurational equilibrium distribution to a distribution that can be easily sampled. Accurate computation of free-energy differences and discovery of new configurations are demonstrated, providing a statistical mechanics tool that can avoid rare events during sampling without prior knowledge of reaction coordinates.
Collapse
|
35
|
Abstract
![]()
A correct estimate
of ligand binding modes and a ratio of their
occupancies is crucial for calculations of binding free energies.
The newly developed method BLUES combines molecular dynamics with
nonequilibrium candidate Monte Carlo. Nonequilibrium candidate Monte
Carlo generates a plethora of possible binding modes and molecular
dynamics enables the system to relax. We used BLUES to investigate
binding modes of caffeine in the active site of its metabolizing enzyme
Cytochrome P450 1A2 with the aim of elucidating metabolite-formation
profiles at different concentrations. Because the activation energies
of all sites of metabolism do not show a clear preference for one
metabolite over the others, the orientations in the active site must
play a key role. In simulations with caffeine located in a spacious
pocket above the I-helix, it points N3 and N1 to the heme iron, whereas
in simulations where caffeine is in close proximity to the heme N7
and C8 are preferably oriented toward the heme iron. We propose a
mechanism where at low caffeine concentrations caffeine binds to the
upper part of the active site, leading to formation of the main metabolite
paraxanthine. On the other hand, at high concentrations two molecules
are located in the active site, forcing one molecule into close proximity
to the heme and yielding metabolites theophylline and trimethyluretic
acid. Our results offer an explanation of previously published experimental
results.
Collapse
Affiliation(s)
- Zuzana Jandova
- Institute of Molecular Modeling and Simulation , University of Natural Resources and Life Sciences, Vienna , 1180 Vienna , Austria
| | | | | | | | - Chris Oostenbrink
- Institute of Molecular Modeling and Simulation , University of Natural Resources and Life Sciences, Vienna , 1180 Vienna , Austria
| |
Collapse
|
36
|
Spiriti J, Subramanian SR, Palli R, Wu M, Zuckerman DM. Middle-way flexible docking: Pose prediction using mixed-resolution Monte Carlo in estrogen receptor α. PLoS One 2019; 14:e0215694. [PMID: 31013302 PMCID: PMC6478315 DOI: 10.1371/journal.pone.0215694] [Citation(s) in RCA: 5] [Impact Index Per Article: 1.0] [Reference Citation Analysis] [Abstract] [MESH Headings] [Grants] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 01/31/2019] [Accepted: 04/06/2019] [Indexed: 12/17/2022] Open
Abstract
There is a vast gulf between the two primary strategies for simulating protein-ligand interactions. Docking methods significantly limit or eliminate protein flexibility to gain great speed at the price of uncontrolled inaccuracy, whereas fully flexible atomistic molecular dynamics simulations are expensive and often suffer from limited sampling. We have developed a flexible docking approach geared especially for highly flexible or poorly resolved targets based on mixed-resolution Monte Carlo (MRMC), which is intended to offer a balance among speed, protein flexibility, and sampling power. The binding region of the protein is treated with a standard atomistic force field, while the remainder of the protein is modeled at the residue level with a Gō model that permits protein flexibility while saving computational cost. Implicit solvation is used. Here we assess three facets of the MRMC approach with implications for other docking studies: (i) the role of receptor flexibility in cross-docking pose prediction; (ii) the use of non-equilibrium candidate Monte Carlo (NCMC) and (iii) the use of pose-clustering in scoring. We examine 61 co-crystallized ligands of estrogen receptor α, an important cancer target known for its flexibility. We also compare the performance of the MRMC approach with Autodock smina. Adding protein flexibility, not surprisingly, leads to significantly lower total energies and stronger interactions between protein and ligand, but notably we document the important role of backbone flexibility in the improvement. The improved backbone flexibility also leads to improved performance relative to smina. Somewhat unexpectedly, our implementation of NCMC leads to only modestly improved sampling of ligand poses. Overall, the addition of protein flexibility improves the performance of docking, as measured by energy-ranked poses, but we do not find significant improvements based on cluster information or the use of NCMC. We discuss possible improvements for the model including alternative coarse-grained force fields, improvements to the treatment of solvation, and adding additional types of NCMC moves.
Collapse
Affiliation(s)
- Justin Spiriti
- Department of Biomedical Engineering, Oregon Health and Science University, Portland, OR 97239, United States of America
| | - Sundar Raman Subramanian
- Department of Computational and Systems Biology, University of Pittsburgh, Pittsburgh, PA 15260, United States of America
| | - Rohith Palli
- Department of Computational and Systems Biology, University of Pittsburgh, Pittsburgh, PA 15260, United States of America
| | - Maria Wu
- Department of Computational and Systems Biology, University of Pittsburgh, Pittsburgh, PA 15260, United States of America
| | - Daniel M. Zuckerman
- Department of Biomedical Engineering, Oregon Health and Science University, Portland, OR 97239, United States of America
- * E-mail:
| |
Collapse
|
37
|
Rotskoff GM, Vanden-Eijnden E. Dynamical Computation of the Density of States and Bayes Factors Using Nonequilibrium Importance Sampling. PHYSICAL REVIEW LETTERS 2019; 122:150602. [PMID: 31050526 DOI: 10.1103/physrevlett.122.150602] [Citation(s) in RCA: 6] [Impact Index Per Article: 1.2] [Reference Citation Analysis] [Abstract] [Track Full Text] [Subscribe] [Scholar Register] [Received: 10/01/2018] [Revised: 02/27/2019] [Indexed: 06/09/2023]
Abstract
Nonequilibrium sampling is potentially much more versatile than its equilibrium counterpart, but it comes with challenges because the invariant distribution is not typically known when the dynamics breaks detailed balance. Here, we derive a generic importance sampling technique that leverages the statistical power of configurations transported by nonequilibrium trajectories and can be used to compute averages with respect to arbitrary target distributions. As a dissipative reweighting scheme, the method can be viewed in relation to the annealed importance sampling (AIS) method and the related Jarzynski equality. Unlike AIS, our approach gives an unbiased estimator, with a provably lower variance than directly estimating the average of an observable. We also establish a direct relation between a dynamical quantity, the dissipation, and the volume of phase space, from which we can compute quantities such as the density of states and Bayes factors. We illustrate the properties of estimators relying on this sampling technique in the context of density of state calculations, showing that it scales favorable with dimensionality-in particular, we show that it can be used to compute the phase diagram of the mean-field Ising model from a single nonequilibrium trajectory. We also demonstrate the robustness and efficiency of the approach with an application to a Bayesian model comparison problem of the type encountered in astrophysics and machine learning.
Collapse
Affiliation(s)
- Grant M Rotskoff
- Courant Institute, New York University, 251 Mercer Street, New York, New York 10012, USA
| | - Eric Vanden-Eijnden
- Courant Institute, New York University, 251 Mercer Street, New York, New York 10012, USA
| |
Collapse
|
38
|
Burley KH, Gill SC, Lim NM, Mobley DL. Enhancing Side Chain Rotamer Sampling Using Nonequilibrium Candidate Monte Carlo. J Chem Theory Comput 2019; 15:1848-1862. [PMID: 30677291 DOI: 10.1021/acs.jctc.8b01018] [Citation(s) in RCA: 15] [Impact Index Per Article: 3.0] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/19/2022]
Abstract
Molecular simulations are a valuable tool for studying biomolecular motions and thermodynamics. However, such motions can be slow compared to simulation time scales, yet critical. Specifically, adequate sampling of side chain motions in protein binding pockets is crucial for obtaining accurate estimates of ligand binding free energies from molecular simulations. The time scale of side chain rotamer flips can range from a few ps to several hundred ns or longer, particularly in crowded environments like the interior of proteins. Here, we apply a mixed nonequilibrium candidate Monte Carlo (NCMC)/molecular dynamics (MD) method to enhance sampling of side chain rotamers. The NCMC portion of our method applies a switching protocol wherein the steric and electrostatic interactions between target side chain atoms and the surrounding environment are cycled off and then back on during the course of a move proposal. Between NCMC move proposals, simulation of the system continues via traditional molecular dynamics. Here, we first validate this approach on a simple, solvated valine-alanine dipeptide system and then apply it to a well-studied model ligand binding site in T4 lysozyme L99A. We compute the rate of rotamer transitions for a valine side chain using our approach and compare it to that of traditional molecular dynamics simulations. Here, we show that our NCMC/MD method substantially enhances side chain sampling, especially in systems where the torsional barrier to rotation is high (≥10 kcal/mol). These barriers can be intrinsic torsional barriers or steric barriers imposed by the environment. Overall, this may provide a promising strategy to selectively improve side chain sampling in molecular simulations.
Collapse
Affiliation(s)
- Kalistyn H Burley
- Department of Pharmaceutical Sciences , University of California, Irvine , Irvine , California 92697 , United States
| | - Samuel C Gill
- Department of Chemistry , University of California, Irvine , Irvine , California 92617 , United States
| | - Nathan M Lim
- Department of Pharmaceutical Sciences , University of California, Irvine , Irvine , California 92697 , United States
| | - David L Mobley
- Department of Pharmaceutical Sciences , University of California, Irvine , Irvine , California 92697 , United States.,Department of Chemistry , University of California, Irvine , Irvine , California 92617 , United States
| |
Collapse
|
39
|
Gilabert JF, Lecina D, Estrada J, Guallar V. Monte Carlo Techniques for Drug Design: The Success Case of PELE. BIOMOLECULAR SIMULATIONS IN STRUCTURE-BASED DRUG DISCOVERY 2018. [DOI: 10.1002/9783527806836.ch5] [Citation(s) in RCA: 8] [Impact Index Per Article: 1.3] [Reference Citation Analysis] [Track Full Text] [Subscribe] [Scholar Register] [Indexed: 12/20/2022]
Affiliation(s)
- Joan F. Gilabert
- Barcelona Supercomputing Center (BSC); Life Science Department; Jordi Girona 29 08034 Barcelona Spain
| | - Daniel Lecina
- Barcelona Supercomputing Center (BSC); Life Science Department; Jordi Girona 29 08034 Barcelona Spain
| | - Jorge Estrada
- Barcelona Supercomputing Center (BSC); Life Science Department; Jordi Girona 29 08034 Barcelona Spain
| | - Victor Guallar
- Barcelona Supercomputing Center (BSC); Life Science Department; Jordi Girona 29 08034 Barcelona Spain
- ICREA; Passeig Lluís Companys 23 08010 Barcelona Spain
| |
Collapse
|
40
|
Rathee VS, Sidky H, Sikora BJ, Whitmer JK. Role of Associative Charging in the Entropy-Energy Balance of Polyelectrolyte Complexes. J Am Chem Soc 2018; 140:15319-15328. [PMID: 30351015 DOI: 10.1021/jacs.8b08649] [Citation(s) in RCA: 65] [Impact Index Per Article: 10.8] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 02/08/2023]
Abstract
Polyelectrolytes may be classified into two primary categories (strong and weak) depending on how their charge state responds to the local environment. Both of these find use in many applications, including drug delivery, gene therapy, layer-by-layer films, and fabrication of ion filtration membranes. The mechanism of polyelectrolyte complexation is, however, still not completely understood, though experimental investigations suggest that entropy gain due to release of counterions is the key driving force for strong polyelectrolyte complexation. Here we perform a comprehensive thermodynamic investigation through coarse-grained molecular simulations permitting us to calculate the free energy of complex formation. Importantly, our expanded-ensemble methods permit the explicit separation of energetic and entropic contributions to the free energy. Our investigations indicate that entropic contributions indeed dominate the free energy of complex formation for strong polyelectrolytes, but are less important than energetic contributions when weak electrostatic coupling or weak polyelectrolytes are present. Our results provide a new view of the free energy of polyelectrolyte complex formation driven by polymer association, which should also arise in systems with large charge spacings or bulky counterions, both of which act to weaken ion-polymer binding.
Collapse
Affiliation(s)
- Vikramjit S Rathee
- Department of Chemical and Biomolecular Engineering , University of Notre Dame , Notre Dame , Indiana 46556 , United States
| | - Hythem Sidky
- Department of Chemical and Biomolecular Engineering , University of Notre Dame , Notre Dame , Indiana 46556 , United States
| | - Benjamin J Sikora
- Department of Chemical and Biomolecular Engineering , University of Notre Dame , Notre Dame , Indiana 46556 , United States
| | - Jonathan K Whitmer
- Department of Chemical and Biomolecular Engineering , University of Notre Dame , Notre Dame , Indiana 46556 , United States
| |
Collapse
|
41
|
Gill SC, Lim NM, Grinaway PB, Rustenburg AS, Fass J, Ross GA, Chodera JD, Mobley DL. Binding Modes of Ligands Using Enhanced Sampling (BLUES): Rapid Decorrelation of Ligand Binding Modes via Nonequilibrium Candidate Monte Carlo. J Phys Chem B 2018; 122:5579-5598. [PMID: 29486559 PMCID: PMC5980761 DOI: 10.1021/acs.jpcb.7b11820] [Citation(s) in RCA: 50] [Impact Index Per Article: 8.3] [Reference Citation Analysis] [Abstract] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/20/2023]
Abstract
Accurately predicting protein-ligand binding affinities and binding modes is a major goal in computational chemistry, but even the prediction of ligand binding modes in proteins poses major challenges. Here, we focus on solving the binding mode prediction problem for rigid fragments. That is, we focus on computing the dominant placement, conformation, and orientations of a relatively rigid, fragment-like ligand in a receptor, and the populations of the multiple binding modes which may be relevant. This problem is important in its own right, but is even more timely given the recent success of alchemical free energy calculations. Alchemical calculations are increasingly used to predict binding free energies of ligands to receptors. However, the accuracy of these calculations is dependent on proper sampling of the relevant ligand binding modes. Unfortunately, ligand binding modes may often be uncertain, hard to predict, and/or slow to interconvert on simulation time scales, so proper sampling with current techniques can require prohibitively long simulations. We need new methods which dramatically improve sampling of ligand binding modes. Here, we develop and apply a nonequilibrium candidate Monte Carlo (NCMC) method to improve sampling of ligand binding modes. In this technique, the ligand is rotated and subsequently allowed to relax in its new position through alchemical perturbation before accepting or rejecting the rotation and relaxation as a nonequilibrium Monte Carlo move. When applied to a T4 lysozyme model binding system, this NCMC method shows over 2 orders of magnitude improvement in binding mode sampling efficiency compared to a brute force molecular dynamics simulation. This is a first step toward applying this methodology to pharmaceutically relevant binding of fragments and, eventually, drug-like molecules. We are making this approach available via our new Binding modes of ligands using enhanced sampling (BLUES) package which is freely available on GitHub.
Collapse
Affiliation(s)
- Samuel C. Gill
- Department of Chemistry, University of California, Irvine
| | - Nathan M. Lim
- Department of Pharmaceutical Sciences, University of California, Irvine
| | - Patrick B. Grinaway
- Graduate Program in Physiology, Biophysics, and Systems Biology, Weill Cornell Medical College, New York, NY 10065
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065
| | - Ariën S. Rustenburg
- Graduate Program in Physiology, Biophysics, and Systems Biology, Weill Cornell Medical College, New York, NY 10065
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065
| | - Josh Fass
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065
- Tri-Institutional PhD Program in Computational Biology and Medicine, New York, NY 10065
| | - Gregory A. Ross
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065
| | | | - David L. Mobley
- Department of Chemistry, University of California, Irvine
- Department of Pharmaceutical Sciences, University of California, Irvine
| |
Collapse
|
42
|
Ross GA, Rustenburg AS, Grinaway PB, Fass J, Chodera JD. Biomolecular Simulations under Realistic Macroscopic Salt Conditions. J Phys Chem B 2018; 122:5466-5486. [PMID: 29649876 PMCID: PMC6078207 DOI: 10.1021/acs.jpcb.7b11734] [Citation(s) in RCA: 62] [Impact Index Per Article: 10.3] [Reference Citation Analysis] [Abstract] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/06/2023]
Abstract
Biomolecular simulations are typically performed in an aqueous environment where the number of ions remains fixed for the duration of the simulation, generally with either a minimally neutralizing ion environment or a number of salt pairs intended to match the macroscopic salt concentration. In contrast, real biomolecules experience local ion environments where the salt concentration is dynamic and may differ from bulk. The degree of salt concentration variability and average deviation from the macroscopic concentration remains, as yet, unknown. Here, we describe the theory and implementation of a Monte Carlo osmostat that can be added to explicit solvent molecular dynamics or Monte Carlo simulations to sample from a semigrand canonical ensemble in which the number of salt pairs fluctuates dynamically during the simulation. The osmostat reproduces the correct equilibrium statistics for a simulation volume that can exchange ions with a large reservoir at a defined macroscopic salt concentration. To achieve useful Monte Carlo acceptance rates, the method makes use of nonequilibrium candidate Monte Carlo (NCMC) moves in which monovalent ions and water molecules are alchemically transmuted using short nonequilibrium trajectories, with a modified Metropolis-Hastings criterion ensuring correct equilibrium statistics for an ( Δμ, N, p, T) ensemble to achieve a ∼1046× boost in acceptance rates. We demonstrate how typical protein (DHFR and the tyrosine kinase Src) and nucleic acid (Drew-Dickerson B-DNA dodecamer) systems exhibit salt concentration distributions that significantly differ from fixed-salt bulk simulations and display fluctuations that are on the same order of magnitude as the average.
Collapse
Affiliation(s)
- Gregory A. Ross
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065
- Present address: Schrödinger, New York, NY 10036
| | - Ariën S. Rustenburg
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065
- Graduate Program in Physiology, Biophysics, and Systems Biology, Weill Cornell Medical College, New York, NY 10065
| | - Patrick B. Grinaway
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065
- Graduate Program in Physiology, Biophysics, and Systems Biology, Weill Cornell Medical College, New York, NY 10065
| | - Josh Fass
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065
- Tri-Institutional Training Program in Computational Biology and Medicine, New York, NY 10065
| | - John D. Chodera
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065
| |
Collapse
|
43
|
Iglesias J, Saen‐oon S, Soliva R, Guallar V. Computational structure‐based drug design: Predicting target flexibility. WILEY INTERDISCIPLINARY REVIEWS-COMPUTATIONAL MOLECULAR SCIENCE 2018. [DOI: 10.1002/wcms.1367] [Citation(s) in RCA: 11] [Impact Index Per Article: 1.8] [Reference Citation Analysis] [Track Full Text] [Subscribe] [Scholar Register] [Indexed: 01/27/2023]
Affiliation(s)
| | | | | | - Victor Guallar
- Life Science DepartmentBarcelonaSpain
- ICREA, Passeig Lluís Companys 23BarcelonaSpain
| |
Collapse
|
44
|
Fass J, Sivak DA, Crooks GE, Beauchamp KA, Leimkuhler B, Chodera JD. Quantifying Configuration-Sampling Error in Langevin Simulations of Complex Molecular Systems. ENTROPY 2018; 20. [PMID: 30393452 PMCID: PMC6208357 DOI: 10.3390/e20050318] [Citation(s) in RCA: 19] [Impact Index Per Article: 3.2] [Reference Citation Analysis] [Abstract] [Key Words] [Grants] [Track Full Text] [Download PDF] [Figures] [Subscribe] [Scholar Register] [Indexed: 11/25/2022]
Abstract
While Langevin integrators are popular in the study of equilibrium properties of complex systems, it is challenging to estimate the timestep-induced discretization error: the degree to which the sampled phase-space or configuration-space probability density departs from the desired target density due to the use of a finite integration timestep. Sivak et al., introduced a convenient approach to approximating a natural measure of error between the sampled density and the target equilibrium density, the Kullback-Leibler (KL) divergence, in phase space, but did not specifically address the issue of configuration-space properties, which are much more commonly of interest in molecular simulations. Here, we introduce a variant of this near-equilibrium estimator capable of measuring the error in the configuration-space marginal density, validating it against a complex but exact nested Monte Carlo estimator to show that it reproduces the KL divergence with high fidelity. To illustrate its utility, we employ this new near-equilibrium estimator to assess a claim that a recently proposed Langevin integrator introduces extremely small configuration-space density errors up to the stability limit at no extra computational expense. Finally, we show how this approach to quantifying sampling bias can be applied to a wide variety of stochastic integrators by following a straightforward procedure to compute the appropriate shadow work, and describe how it can be extended to quantify the error in arbitrary marginal or conditional distributions of interest.
Collapse
Affiliation(s)
- Josh Fass
- Tri-Institutional PhD Program in Computational Biology & Medicine, New York, NY 10065, USA
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065, USA
| | - David A. Sivak
- Department of Physics, Simon Fraser University, Burnaby, BC V5A 1S6, Canada
| | | | | | - Benedict Leimkuhler
- School of Mathematics and Maxwell Institute of Mathematical Sciences, James Clerk Maxwell Building, Kings Buildings, University of Edinburgh, Edinburgh EH9 3FD, UK
| | - John D. Chodera
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065, USA
- Correspondence:
| |
Collapse
|
45
|
Suh D, Radak BK, Chipot C, Roux B. Enhanced configurational sampling with hybrid non-equilibrium molecular dynamics–Monte Carlo propagator. J Chem Phys 2018; 148:014101. [DOI: 10.1063/1.5004154] [Citation(s) in RCA: 17] [Impact Index Per Article: 2.8] [Reference Citation Analysis] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/10/2023] Open
Affiliation(s)
- Donghyuk Suh
- Department of Chemistry, University of Chicago, Chicago, Illinois 60637, USA
| | - Brian K. Radak
- Leadership Computing Facility, Argonne National Laboratory, Argonne, Illinois 60439-8643, USA
| | - Christophe Chipot
- Laboratoire International Associé Centre National de la Recherche Scientifique and University of Illinois at Urbana-Champaign, UMR 7565, Université de Lorraine, BP 70239, 54506 Vandævre-lès-Nancy, France
- Department of Physics, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA and Beckman Institute, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA
| | - Benoît Roux
- Department of Biochemistry and Molecular Biology, University of Chicago, Chicago, Illinois 60637-1454, USA
- Center for Nanoscale Materials, Argonne National Laboratory, Argonne, Illinois 60439-8643, USA
| |
Collapse
|
46
|
Kurut A, Fonseca R, Boomsma W. Driving Structural Transitions in Molecular Simulations Using the Nonequilibrium Candidate Monte Carlo. J Phys Chem B 2018; 122:1195-1204. [PMID: 29260565 DOI: 10.1021/acs.jpcb.7b11426] [Citation(s) in RCA: 5] [Impact Index Per Article: 0.8] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 11/29/2022]
Abstract
Hybrid simulation procedures which combine molecular dynamics with Monte Carlo are attracting increasing attention as tools for improving the sampling efficiency in molecular simulations. In particular, encouraging results have been reported for nonequilibrium candidate protocols, in which a Monte Carlo move is applied gradually, and interleaved with a process that equilibrates the remaining degrees of freedom. Although initial studies have uncovered a substantial potential of the method, its practical applicability for sampling structural transitions in macromolecules remains incompletely understood. Here, we address this issue by systematically investigating the efficiency of the nonequilibrium candidate Monte Carlo on the sampling of rotameric distributions of two peptide systems at atomistic resolution both in vacuum and explicit solvent. The studied systems allow us to directly probe the efficiency with which a single or a few slow degrees of freedom can be driven between well-separated free-energy minima and to explore the sensitivity of the method toward the involved free parameters. In line with results on other systems, our study suggests that order-of-magnitude gains can be obtained in certain scenarios but also identifies challenges that arise when applying the procedure in explicit solvent.
Collapse
Affiliation(s)
- Anıl Kurut
- Department of Computer Science, University of Copenhagen , 2100 Copenhagen Ø, Denmark
| | - Rasmus Fonseca
- Department of Molecular and Cellular Physiology, Stanford University , Stanford, California 94305, United States
| | - Wouter Boomsma
- Department of Computer Science, University of Copenhagen , 2100 Copenhagen Ø, Denmark
| |
Collapse
|
47
|
Radak BK, Chipot C, Suh D, Jo S, Jiang W, Phillips JC, Schulten K, Roux B. Constant-pH Molecular Dynamics Simulations for Large Biomolecular Systems. J Chem Theory Comput 2017; 13:5933-5944. [PMID: 29111720 DOI: 10.1021/acs.jctc.7b00875] [Citation(s) in RCA: 123] [Impact Index Per Article: 17.6] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 02/04/2023]
Abstract
An increasingly important endeavor is to develop computational strategies that enable molecular dynamics (MD) simulations of biomolecular systems with spontaneous changes in protonation states under conditions of constant pH. The present work describes our efforts to implement the powerful constant-pH MD simulation method, based on a hybrid nonequilibrium MD/Monte Carlo (neMD/MC) technique within the highly scalable program NAMD. The constant-pH hybrid neMD/MC method has several appealing features; it samples the correct semigrand canonical ensemble rigorously, the computational cost increases linearly with the number of titratable sites, and it is applicable to explicit solvent simulations. The present implementation of the constant-pH hybrid neMD/MC in NAMD is designed to handle a wide range of biomolecular systems with no constraints on the choice of force field. Furthermore, the sampling efficiency can be adaptively improved on-the-fly by adjusting algorithmic parameters during the simulation. Illustrative examples emphasizing medium- and large-scale applications on next-generation supercomputing architectures are provided.
Collapse
Affiliation(s)
- Brian K Radak
- Leadership Computing Facility, Argonne National Laboratory , Argonne, Illinois 60439-8643, United States
| | - Christophe Chipot
- Laboratoire International Associé Centre National de la Recherche Scientifique et University of Illinois at Urbana-Champaign, Unité Mixte de Recherche No. 7565, Université de Lorraine, Université de Lorraine , B.P. 70239, 54506 Vandœuvre-lès-Nancy cedex, France.,Department of Physics, University of Illinois at Urbana-Champaign , Urbana, Illinois 61801-2325, United States
| | - Donghyuk Suh
- Department of Chemistry, University of Chicago , Chicago, Illinois 60637-1454, United States
| | - Sunhwan Jo
- Leadership Computing Facility, Argonne National Laboratory , Argonne, Illinois 60439-8643, United States
| | - Wei Jiang
- Leadership Computing Facility, Argonne National Laboratory , Argonne, Illinois 60439-8643, United States
| | - James C Phillips
- Theoretical and Computational Biophysics Group, Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign , Urbana, Illinois 61801-2325, United States
| | - Klaus Schulten
- Department of Physics, University of Illinois at Urbana-Champaign , Urbana, Illinois 61801-2325, United States.,Theoretical and Computational Biophysics Group, Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign , Urbana, Illinois 61801-2325, United States
| | - Benoît Roux
- Department of Chemistry, University of Chicago , Chicago, Illinois 60637-1454, United States.,Department of Biochemistry and Molecular Biology, University of Chicago , Chicago, Illinois 60637-1454, United States.,Center for Nanoscale Materials, Argonne National Laboratory , Argonne, Illinois 60439-8643, United States
| |
Collapse
|
48
|
Abstract
An important limitation of standard classical molecular dynamics simulations is the inability to make or break chemical bonds. This restricts severely our ability to study processes that involve even the simplest of chemical reactions, the transfer of a proton. Existing approaches for allowing proton transfer in the context of classical mechanics are rather cumbersome and have not achieved widespread use and routine status. Here we reconsider the combination of molecular dynamics with periodic stochastic proton hops. To ensure computational efficiency, we propose a non-Boltzmann acceptance criterion that is heuristically adjusted to maintain the correct or desirable thermodynamic equilibria between different protonation states and proton transfer rates. Parameters are proposed for hydronium, Asp, Glu, and His. The algorithm is implemented in the program CHARMM and tested on proton diffusion in bulk water and carbon nanotubes and on proton conductance in the gramicidin A channel. Using hopping parameters determined from proton diffusion in bulk water, the model reproduces the enhanced proton diffusivity in carbon nanotubes and gives a reasonable estimate of the proton conductance in gramicidin A.
Collapse
Affiliation(s)
- Themis Lazaridis
- Department of Chemistry, City College of New York/CUNY , 160 Convent Avenue, New York, New York 10031, United States.,Graduate Programs in Chemistry, Biochemistry & Physics, Graduate Center, City University of New York , 365 Fifth Ave, New York, New York 10016, United States
| | - Gerhard Hummer
- Department of Theoretical Biophysics, Max Planck Institute of Biophysics , Max-von-Laue Strasse 3, 60438 Frankfurt am Main, Germany
| |
Collapse
|
49
|
Eastman P, Swails J, Chodera JD, McGibbon RT, Zhao Y, Beauchamp KA, Wang LP, Simmonett AC, Harrigan MP, Stern CD, Wiewiora RP, Brooks BR, Pande VS. OpenMM 7: Rapid development of high performance algorithms for molecular dynamics. PLoS Comput Biol 2017; 13:e1005659. [PMID: 28746339 PMCID: PMC5549999 DOI: 10.1371/journal.pcbi.1005659] [Citation(s) in RCA: 1373] [Impact Index Per Article: 196.1] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 10/20/2016] [Revised: 08/09/2017] [Accepted: 06/27/2017] [Indexed: 01/22/2023] Open
Abstract
OpenMM is a molecular dynamics simulation toolkit with a unique focus on extensibility. It allows users to easily add new features, including forces with novel functional forms, new integration algorithms, and new simulation protocols. Those features automatically work on all supported hardware types (including both CPUs and GPUs) and perform well on all of them. In many cases they require minimal coding, just a mathematical description of the desired function. They also require no modification to OpenMM itself and can be distributed independently of OpenMM. This makes it an ideal tool for researchers developing new simulation methods, and also allows those new methods to be immediately available to the larger community.
Collapse
Affiliation(s)
- Peter Eastman
- Department of Chemistry, Stanford University, Stanford, California, United States of America
| | - Jason Swails
- Department of Chemistry and Chemical Biology and BioMaPS Institute, Rutgers University, Piscataway, New Jersey, United States of America
| | - John D Chodera
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, New York, United States of America
| | - Robert T McGibbon
- Department of Chemistry, Stanford University, Stanford, California, United States of America
| | - Yutong Zhao
- Department of Chemistry, Stanford University, Stanford, California, United States of America
| | - Kyle A Beauchamp
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, New York, United States of America
| | - Lee-Ping Wang
- Department of Chemistry, University of California, Davis, Davis, California, United States of America
| | - Andrew C Simmonett
- Laboratory of Computational Biology, National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, Maryland, United States of America
| | - Matthew P Harrigan
- Department of Chemistry, Stanford University, Stanford, California, United States of America
| | - Chaya D Stern
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, New York, United States of America.,Tri-Institutional PhD Program in Chemical Biology, Memorial Sloan Kettering Cancer Center, New York, New York, United States of America
| | - Rafal P Wiewiora
- Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, New York, United States of America.,Tri-Institutional PhD Program in Chemical Biology, Memorial Sloan Kettering Cancer Center, New York, New York, United States of America
| | - Bernard R Brooks
- Laboratory of Computational Biology, National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, Maryland, United States of America
| | - Vijay S Pande
- Department of Chemistry, Stanford University, Stanford, California, United States of America.,Department of Computer Science, Stanford University, Stanford, California, United States of America
| |
Collapse
|
50
|
Radak BK, Roux B. Efficiency in nonequilibrium molecular dynamics Monte Carlo simulations. J Chem Phys 2017; 145:134109. [PMID: 27782441 DOI: 10.1063/1.4964288] [Citation(s) in RCA: 15] [Impact Index Per Article: 2.1] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/11/2022] Open
Abstract
Hybrid algorithms combining nonequilibrium molecular dynamics and Monte Carlo (neMD/MC) offer a powerful avenue for improving the sampling efficiency of computer simulations of complex systems. These neMD/MC algorithms are also increasingly finding use in applications where conventional approaches are impractical, such as constant-pH simulations with explicit solvent. However, selecting an optimal nonequilibrium protocol for maximum efficiency often represents a non-trivial challenge. This work evaluates the efficiency of a broad class of neMD/MC algorithms and protocols within the theoretical framework of linear response theory. The approximations are validated against constant pH-MD simulations and shown to provide accurate predictions of neMD/MC performance. An assessment of a large set of protocols confirms (both theoretically and empirically) that a linear work protocol gives the best neMD/MC performance. Finally, a well-defined criterion for optimizing the time parameters of the protocol is proposed and demonstrated with an adaptive algorithm that improves the performance on-the-fly with minimal cost.
Collapse
Affiliation(s)
- Brian K Radak
- Leadership Computing Facility, Argonne National Laboratory, Argonne, Illinois 60439-8643, USA
| | - Benoît Roux
- Department of Biochemistry and Molecular Biology, University of Chicago, Chicago, Illinois 60637-1454, USA
| |
Collapse
|