Temperature-accelerated Molecular Dynamics Research Articles

Formation of solid electrolyte interphase in Li-ion batteries: Insights from temperature-accelerated ReaxFF molecular dynamics

The solid electrolyte interphase (SEI) forms between anode and electrolyte in the lithium-ion battery (LIB) and is considered to be the main contributor to LIB degradation. To understand the degradation and SEI formation, we use reactive force-field (ReaxFF) molecular dynamics to gain atomic level insights into the key phenomena. We report the effect of different reaction temperatures on SEI formation on the lithium metal anode surface in contact with pure ethylene carbonate electrolyte. The temperature elevation is found to boost the SEI formation. To investigate the structural variations of SEI, we compare the quantitative variations of key species and perform analysis using the charge distribution and radial distribution function (RDF). Our work reveals that elevated temperature of up to 700 K is appropriate for acceleration of chemical reactions during SEI formation. The study provides insights into structural variations of the lithium metal electrode and SEI at the molecular level, demonstrating that temperature-accelerated reactive molecular dynamics is an efficient and effective tool for predicting the performance and degradation of the Li-ion batteries.

Unbiasing Enhanced Sampling on a High-Dimensional Free Energy Surface with a Deep Generative Model.

Biased enhanced sampling methods that utilize collective variables (CVs) are powerful tools for sampling conformational ensembles. Due to their large intrinsic dimensions, efficiently generating conformational ensembles for complex systems requires enhanced sampling on high-dimensional free energy surfaces. While temperature-accelerated molecular dynamics (TAMD) can trivially adopt many CVs in a simulation, unbiasing the simulation to generate unbiased conformational ensembles requires accurate modeling of a high-dimensional CV probability distribution, which is challenging for traditional density estimation techniques. Here we propose an unbiasing method based on the score-based diffusion model, a deep generative learning method that excels in density estimation across complex data landscapes. We demonstrate that this unbiasing approach, tested on multiple TAMD simulations, significantly outperforms traditional unbiasing methods and can generate accurate unbiased conformational ensembles. With the proposed approach, TAMD can adopt CVs that focus on improving sampling efficiency and the proposed unbiasing method enables accurate evaluation of ensemble averages of important chemical features.

An interoperable implementation of collective-variable based enhanced sampling methods in extended phase space within the OpenMM package.

Collective variable (CV)-based enhanced sampling techniques are widely used today for accelerating barrier-crossing events in molecular simulations. A class of these methods, which includes temperature accelerated molecular dynamics (TAMD)/driven-adiabatic free energy dynamics (d-AFED), unified free energy dynamics (UFED), and temperature accelerated sliced sampling (TASS), uses an extended variable formalism to achieve quick exploration of conformational space. These techniques are powerful, as they enhance the sampling of a large number of CVs simultaneously compared to other techniques. Extended variables are kept at a much higher temperature than the physical temperature by ensuring adiabatic separation between the extended and physical subsystems and employing rigorous thermostatting. In this work, we present a computational platform to perform extended phase space enhanced sampling simulations using the open-source molecular dynamics engine OpenMM. The implementation allows users to have interoperability of sampling techniques, as well as employ state-of-the-art thermostats and multiple time-stepping. This work also presents protocols for determining the critical parameters and procedures for reconstructing high-dimensional free energy surfaces. As a demonstration, we present simulation results on the high dimensional conformational landscapes of the alanine tripeptide in vacuo, tetra-N-methylglycine (tetra-sarcosine) peptoid in implicit solvent, and the Trp-cage mini protein in explicit water.

Molecular Dynamics Simulations of Ion Permeation in Human Voltage-Gated Sodium Channels.

The recent determination of cryo-EM structures of voltage-gated sodium (Nav) channels has revealed many details of these proteins. However, knowledge of ionic permeation through the Nav pore remains limited. In this work, we performed atomistic molecular dynamics (MD) simulations to study the structural features of various neuronal Nav channels based on homology modeling of the cryo-EM structure of the human Nav1.4 channel and, in addition, on the recently resolved configuration for Nav1.2. In particular, single Na+ permeation events during standard MD runs suggest that the ion resides in the inner part of the Nav selectivity filter (SF). On-the-fly free energy parametrization (OTFP) temperature-accelerated molecular dynamics (TAMD) was also used to calculate two-dimensional free energy surfaces (FESs) related to single/double Na+ translocation through the SF of the homology-based Nav1.2 model and the cryo-EM Nav1.2 structure, with different realizations of the DEKA filter domain. These additional simulations revealed distinct mechanisms for single and double Na+ permeation through the wild-type SF, which has a charged lysine in the DEKA ring. Moreover, the configurations of the ions in the SF corresponding to the metastable states of the FESs are specific for each SF motif. Overall, the description of these mechanisms gives us new insights into ion conduction in human Nav cryo-EM-based and cryo-EM configurations that could advance understanding of these systems and how they differ from potassium and bacterial Nav channels.

Improved Sampling and Free Energy Estimates for Antibiotic Permeation through Bacterial Porins.

Antibiotics enter into bacterial cells via protein channels that serve as low-energy pathways through the outer membrane, which is otherwise impenetrable. Insights into the molecular mechanisms underlying the transport processes are vital for the development of effective antibacterials. A much-desired prerequisite is an accurate and reproducible determination of free energy surfaces for antibiotic translocation, enabling quantitative and meaningful comparisons of permeation mechanisms for different classes of antibiotics. Inefficient sampling along the orthogonal degrees of freedom, for example, in umbrella sampling and metadynamics approaches, is however a key limitation affecting the accuracy and the convergence of free energy estimates. To overcome this limitation, two sampling methods have been employed in the present study that, respectively, combine umbrella sampling and metadynamics-style biasing schemes with temperature acceleration for improved sampling along orthogonal degrees of freedom. As a model for the transport of bulky solutes, the ciprofloxacin-OmpF system has been selected. The well-tempered metadynamics approach with multiple walkers is compared with its "temperature-accelerated" variant in terms of improvements in sampling and convergence of free energy estimates. We find that the inclusion of collective variables governing solute degrees of freedom and solute-water interactions within the sampling scheme largely alleviates sampling issues. Concerning improved sampling and convergence of free energy estimates from independent simulations, the temperature-accelerated sliced sampling approach that combines umbrella sampling with temperature-accelerated molecular dynamics performs even better as shown for the ciprofloxacin-OmpF system.

Analysis of Charged Peptide Loop-Flipping across a Lipid Bilayer Using the String Method with Swarms of Trajectories.

The hydrophobic core of the lipid bilayer is conventionally considered a barrier to the translocation of charged species such that the translocation of even single ions occurs on long time scales. In contrast, experiments have revealed that some materials, including peptides, proteins, and nanoparticles, can translocate multiple charged moieties across the bilayer on experimentally relevant time scales. Understanding the molecular mechanisms underlying this behavior is challenging because resolving corresponding free-energy landscapes with molecular simulation techniques is computationally expensive. To address this challenge, we use atomistic molecular dynamics simulations with the swarms-of-trajectories (SOT) string method to analyze charge translocation pathways across single-component lipid bilayers as a function of multiple collective variables. We first demonstrate that the SOT string method can reproduce the free-energy barrier for the translocation of a charged lysine amino acid analogue in good agreement with the literature. We then obtain minimum free-energy pathways for the translocation, or flipping, of charged peptide loops across the lipid bilayer by utilizing trajectories from prior temperature-accelerated molecular dynamics (TAMD) simulations as initial configurations. The corresponding potential of mean force calculations reveal that the protonation of a central membrane-exposed aspartate residue substantially reduces the free-energy barrier for flipping charged loops by modulating the water content of the bilayer. These results provide new insight into the thermodynamics underlying loop-flipping processes and highlight how the combination of TAMD and the SOT string method can be used to understand complex charge translocation mechanisms.

Characterizing the Translocation of Charged Peptide Loops across Lipid Bilayers with Molecular Dynamics Simulations

The hydrophobic core of the cell membrane is considered largely impermeable to charged molecules because of the large free energy barrier and corresponding long timescale (∼hours) for charge translocation. Contradicting this view, a variety of water-soluble peptides are known to translocate charged groups across the cell membrane on a surprisingly rapid timescale, on the order of few seconds at least for some peptides. In this work, we study the interconversion of a peptide with charged flanking loops between a surface-adsorbed and membrane-embedded state, which requires the translocation (or “flipping”) of charged loops across the bilayer. We utilize all-atom Temperature Accelerated Molecular Dynamics (TAMD) simulations to predict the likelihood of loop flipping without predefining a specific translocation pathway. We show that membrane-exposed charged residues accelerate the flipping of charged peptide loops by stabilizing intramembrane water defects. We further demonstrate that this computational approach can identify multiple flipping pathways without specifying them a priori. The position of charged residues in the transmembrane helix also affects the flipping pathways and the apparent flipping free energy barrier. These detailed molecular-level insights of peptide translocation pathways may be valuable for designing small cationic peptides for applications in gene and drug delivery. Moreover, the methodology and findings discussed here can generalize to more complex behaviors, such as the large-scale conformational rearrangements of integral membrane proteins.

Free energy calculation using space filled design and weighted reconstruction: a modified single sweep approach

ABSTRACTA modified single sweep approach is proposed for generating free energy landscapes. The approach replaces the use of temperature-accelerated molecular dynamics (TAMD) to generate centres in collective variable (CV) space at which mean forces are computed using restrained molecular dynamics (MD) simulations with a sequential space-filling design. This approach also modifies the radial basis function reconstruction step of the traditional single sweep approach and proposes a weighted reconstruction of the free energy surface using the previously generated mean forces. The modified approach is compared to the traditional single sweep (SS) approach on the (φ, ψ) dihedral free-energy map of solvated alanine dipeptide (AD). It is found that the new approach results in a more accurate reconstructed free energy than does the traditional approach when compared to the directly-computed reference free energy landscape. It is shown that the increased accuracy of the overall map stems from the improved 1-dimensional space filling (projective) property of the proposed design compared to that of the TAMD generated centres. A further enhancement in the accuracy of the crucial lower energy regions is enabled by the introduction of weights in the reconstruction step that give more importance to lower energy-valued regions.

Conformational sampling of CpxA: Connecting HAMP motions to the histidine kinase function.

In the histidine kinase family, the HAMP and DHp domains are considered to play an important role into the transmission of signal arising from environmental conditions to the auto-phosphorylation site and to the binding site of response regulator. Several conformational motions inside HAMP have been proposed to transmit this signal: (i) the gearbox model, (ii) α helices rotations, pistons and scissoring, (iii) transition between ordered and disordered states. In the present work, we explore by temperature-accelerated molecular dynamics (TAMD), an enhanced sampling technique, the conformational space of the cytoplasmic region of histidine kinase CpxA. Several HAMP motions, corresponding to α helices rotations, pistoning and scissoring have been detected and correlated to the segmental motions of HAMP and DHp domains of CpxA.

Entry from the Lipid Bilayer: A Possible Pathway for Inhibition of a Peptide G Protein-Coupled Receptor by a Lipophilic Small Molecule.

The pathways that G protein-coupled receptor (GPCR) ligands follow as they bind to or dissociate from their receptors are largely unknown. Protease-activated receptor-1 (PAR1) is a GPCR activated by intramolecular binding of a tethered agonist peptide that is exposed by thrombin cleavage. By contrast, the PAR1 antagonist vorapaxar is a lipophilic drug that binds in a pocket almost entirely occluded from the extracellular solvent. The binding and dissociation pathway of vorapaxar is unknown. Starting with the crystal structure of vorapaxar bound to PAR1, we performed temperature-accelerated molecular dynamics simulations of ligand dissociation. In the majority of simulations, vorapaxar exited the receptor laterally into the lipid bilayer through openings in the transmembrane helix (TM) bundle. Prior to full dissociation, vorapaxar paused in metastable intermediates stabilized by interactions with the receptor and lipid headgroups. Derivatives of vorapaxar with alkyl chains predicted to extend between TM6 and TM7 into the lipid bilayer inhibited PAR1 with apparent on rates similar to that of the parent compound in cell signaling assays. These data are consistent with vorapaxar binding to PAR1 via a pathway that passes between TM6 and TM7 from the lipid bilayer, in agreement with the most consistent pathway observed by molecular dynamics. While there is some evidence of entry of the ligand into rhodopsin and lipid-activated GPCRs from the cell membrane, our study provides the first such evidence for a peptide-activated GPCR and suggests that metastable intermediates along drug binding and dissociation pathways can be stabilized by specific interactions between lipids and the ligand.

Longtime convergence of the temperature-accelerated molecular dynamics method

The equations of the temperature-accelerated molecular dynamics (TAMD) method for the calculations of free energies and partition functions are analyzed. Specifically, the exponential convergence of the law of these stochastic processes is established, with a convergence rate close to the one of the limiting, effective dynamics at higher temperature obtained with infinite acceleration. It is also shown that the invariant measures of TAMD are close to a known reference measure, with an error that can be quantified precisely. Finally, a central limit theorem is proven, which allows the estimation of errors on properties calculated by ergodic time averages. These results not only demonstrate the usefulness and validity range of the TAMD equations, but they also permit in principle to adjust the parameter in these equations to optimize their efficiency.

ELF: An Extended-Lagrangian Free Energy Calculation Module for Multiple Molecular Dynamics Engines.

Extended adaptive biasing force (eABF), a collective variable (CV)-based importance-sampling algorithm, has proven to be very robust and efficient compared with the original ABF algorithm. Its implementation in Colvars, a software addition to molecular dynamics (MD) engines, is, however, currently limited to NAMD and LAMMPS. To broaden the scope of eABF and its variants, like its generalized form (egABF), and make them available to other MD engines, e.g., GROMACS, AMBER, CP2K, and openMM, we present a PLUMED-based implementation, called extended-Lagrangian free energy calculation (ELF). This implementation can be used as a stand-alone gradient estimator for other CV-based sampling algorithms, such as temperature-accelerated MD (TAMD) and extended-Lagrangian metadynamics (MtD). ELF provides the end user with a convenient framework to help select the best-suited importance-sampling algorithm for a given application without any commitment to a particular MD engine.

Self-optimized construction of transition rate matrices from accelerated atomistic simulations with Bayesian uncertainty quantification

A massively parallel method to build large transition rate matrices from temperature-accelerated molecular dynamics trajectories is presented. Bayesian Markov model analysis is used to estimate the expected residence time in the known state space, providing crucial uncertainty quantification for higher-scale simulation schemes such as kinetic Monte Carlo or cluster dynamics. The estimators are additionally used to optimize where exploration is performed and the degree of temperature acceleration on the fly, giving an autonomous, optimal procedure to explore the state space of complex systems. The method is tested against exactly solvable models and used to explore the dynamics of C15 interstitial defects in iron. Our uncertainty quantification scheme allows for accurate modeling of the evolution of these defects over timescales of several seconds.

Effect of Intercalated Water on Potassium Ion Transport through Kv1.2 Channels Studied via On-the-Fly Free-Energy Parametrization.

We introduce a two-dimensional version of the method called on-the-fly free energy parametrization (OTFP) to reconstruct free-energy surfaces using Molecular Dynamics simulations, which we name OTFP-2D. We first test the new method by reconstructing the well-known dihedral angles free energy surface of solvated alanine dipeptide. Then, we use it to investigate the process of K+ ions translocation inside the Kv1.2 channel. By comparing a series of two-dimensional free energy surfaces for ion movement calculated with different conditions on the intercalated water molecules, we first recapitulate the widely accepted knock-on mechanism for ion translocation and then confirm that permeation occurs with water molecules alternated among the ions, in accordance with the latest experimental findings. From a methodological standpoint, our new OTFP-2D algorithm demonstrates the excellent sampling acceleration of temperature-accelerated molecular dynamics and the ability to efficiently compute 2D free-energy surfaces. It will therefore be useful in large variety complex biomacromolecular simulations.

Insight into the Role of the Hv1 C-Terminal Domain in Dimer Stabilization.

The voltage-gated proton-selective channel (Hv1) conducts protons in response to changes in membrane potential. The Hv1 protein forms dimers in the membrane. Crystal structures of Hv1 channels have revealed that the primary contacts between the two monomers are in the C-terminal domain (CTD), which forms a coiled-coil structure. The role of Hv1-CTD in channel assembly and activity is not fully understood. Here, molecular dynamics (MD) simulations of full-length and truncated CTD models of human and mouse Hv1 channels reveal a strong contribution of the CTD to the packing of the transmembrane domains. Simulations of the CTD models highlight four fundamental interactions of the key residues contributing to dimer stability. These include salt bridges, hydrophobic interactions, hydrogen bonds, and a disulfide bond across the dimer interface. At neutral pH, salt-bridge interactions increase dimer stability and the dimer becomes less stable at acidic pH. Hydrophobic core packing of the heptad pattern is important for stability, as shown by favorable nonpolar binding free energies rather than by electrostatic components. Moreover, free-energy calculations indicate that a more uniform hydrophobic core in the coiled-coil structure of the Hv1-NIN, a channel carrying the triple mutation M234N-N235I-V236N, leads to an increase in dimer stability with respect to the wild-type. A Cys disulfide bond has a strong impact on dimer stability by holding the dimer together and facilitating the interactions described above. These results are consistent with dissociative temperatures and energy barriers of dimer dissociation obtained from the temperature-accelerated MD.

Exploring Reaction Pathways for Peptidylprolyl-Isomerase

Pin1 enzyme is one of the peptidyl-prolyl cis/trans isomerases (PPIase), which isomerizes an omega bond of specific phospho-Serine/Threonine-Proline motifs, leading to many biological consequences. However, the isomerization mechanism has not been fully clarified yet though there are several experimental and numerical studies. For example, Tate and coworkers investigated a mutated (C113D) Pin1, and found that the isomerization rate is lower than that in the wild type [1], but it is difficult to understand the mechanism only from the experiment. Previously, Hamelberg and coworkers carried out molecular dynamics (MD) simulations of Pin1 with a model substrate, and because the isomerization process is a rare event, they used accelerated MD method [2] to enhance the sampling in configuration space. We here use temperature accelerated MD (TAMD) with several order parameters [3], and explore the free energy landscape for the isomerization process of Pin1. Combining with the string method [4], we will clarify the reaction pathways, and also examine the reaction rate using milestoning [5] or non-Markov type analysis of trajectories [6].[1] Ning Xu et al. Biochemistry, 53, 5568-5578 (2014).[2] H.A. Velazquez and D. Hamelberg, J. Phys. Chem. B 117, 11509-11517 (2013).[3] H. Fujisaki, K. Moritsugu, Y. Matsunaga, T. Morishita, L. Maragliano, Front. Bioeng. and Biotechnol. 3: 125, Doi: 10.3389/fbioe.2015.00125.[4] L. Maragliano and E. Vanden-Eijnden, Chem. Phys. Lett. 446, 182-190 (2007).[5] L. Maragliano, E. Vanden-Eijnden and B. Roux, J. Chem. Theory Comp. 5, 2589-2594 (2009).[6] E. Suarez, J.L. Adelman, and D.M. Zuckerman, J. Chem. Theor. Comput. DOI:10.1021/acs.jctc.6b00339.

Efficient free energy calculations by combining two complementary tempering sampling methods.

Although energy barriers can be efficiently crossed in the reaction coordinate (RC) guided sampling, this type of method suffers from identification of the correct RCs or requirements of high dimensionality of the defined RCs for a given system. If only the approximate RCs with significant barriers are used in the simulations, hidden energy barriers with small to medium height would exist in other degrees of freedom (DOFs) relevant to the target process and consequently cause the problem of insufficient sampling. To address the sampling in this so-called hidden barrier situation, here we propose an effective approach to combine temperature accelerated molecular dynamics (TAMD), an efficient RC-guided sampling method, with the integrated tempering sampling (ITS), a generalized ensemble sampling method. In this combined ITS-TAMD method, the sampling along the major RCs with high energy barriers is guided by TAMD and the sampling of the rest of the DOFs with lower but not negligible barriers is enhanced by ITS. The performance of ITS-TAMD to three systems in the processes with hidden barriers has been examined. In comparison to the standalone TAMD or ITS approach, the present hybrid method shows three main improvements. (1) Sampling efficiency can be improved at least five times even if in the presence of hidden energy barriers. (2) The canonical distribution can be more accurately recovered, from which the thermodynamic properties along other collective variables can be computed correctly. (3) The robustness of the selection of major RCs suggests that the dimensionality of necessary RCs can be reduced. Our work shows more potential applications of the ITS-TAMD method as the efficient and powerful tool for the investigation of a broad range of interesting cases.

Building Graphs To Describe Dynamics, Kinetics, and Energetics in the d-ALa:d-Lac Ligase VanA.

The d-Ala:d-Lac ligase, VanA, plays a critical role in the resistance of vancomycin. Indeed, it is involved in the synthesis of a peptidoglycan precursor, to which vancomycin cannot bind. The reaction catalyzed by VanA requires the opening of the so-called “ω-loop”, so that the substrates can enter the active site. Here, the conformational landscape of VanA is explored by an enhanced sampling approach: the temperature-accelerated molecular dynamics (TAMD). Analysis of the molecular dynamics (MD) and TAMD trajectories recorded on VanA permits a graphical description of the structural and kinetics aspects of the conformational space of VanA, where the internal mobility and various opening modes of the ω-loop play a major role. The other important feature is the correlation of the ω-loop motion with the movements of the opposite domain, defined as containing the residues A149–Q208. Conformational and kinetic clusters have been determined and a path describing the ω-loop opening was extracted from these clusters. The determination of this opening path, as well as the relative importance of hydrogen bonds along the path, permit one to propose some key residue interactions for the kinetics of the ω-loop opening.

Temperature-accelerated molecular dynamics simulation of the evolution of a low-energy incident Cu3 cluster on the Cu(100) surface with a monoatomic step

The molecular dynamics simulation of the normal incidence of a Cu3 cluster with an energy of 0.2–1 eV/atom on a Cu(100) substrate at an equilibrium temperature of 500–700 K is performed. This substrate contains a rectilinear step, whose height is one atomic layer. After 20-ps relaxation of the atomic cluster on the substrate, its further thermally activated motion on the surface is simulated using the method of temperature- accelerated molecular dynamics. The Cu3 cluster mainly demonstrates rotational-translational motion. The time interval between two successive atomic transitions of the cluster is found to be reduced as the distance to the step is decreased.

Building Graphs to Describe Dynamics, Kinetics and Energetics in the D-Ala:D-Lac Ligase Vana

The D-Ala:D-Lac ligase VanA, plays a very important role into the resistance to vancomycin, as it allows the bacteria to replace the synthesis of D-Ala:D-Ala by the synthesis of D-Ala-D-Lac in the bacterial wall, thus making obsolete the role of vancomycin. The reaction catalyzed by VanA requires the opening of the omega-loop to let the substrates entering the active site. Previous molecular dynamics simulations have shown the spontaneous opening of omega-loop in the presence of a crystallographic disulphide bridge. Here, the conformational landscape of VanA is further explored in the absence of the bridge and by using the enhanced sampling approach temperature-accelerated molecular dynamics (TAMD). This analysis has permitted to reveal a description of the conformational space of VanA, where the internal mobility and various opening modes of omega-loop play a major role. The other important feature is the correlation of the omega-loop motion with the movements of the opposite domain. The relative mobility of these two regions is observed in the conformational space analysis and their coupling is demonstrated by the communities analysis of the protein structure. The observations made on the VanA dynamics and conformations permit to propose some residues interactions to be important for the kinetics of the omega-loop opening. Also, the various graph network analyses performed here provide a panel of methods to describe a protein function in the frame of a larger system-oriented network.