Amber Molecular Dynamics Package Research Articles

Multiscale Methods Framework with the 3D-RISM-KH Molecular Solvation Theory for Supramolecular Structures, Nanomaterials, and Biomolecules: Where Are We Going?

3D-RISM-KH molecular solvation theory based on statistical mechanics has been an engine of the multiscale methods framework, which also includes molecular simulation techniques. Its applications range from the solvation energy of small molecules to the phase behavior of polymers and biomolecules. Molecular solvation theory predicts and explains the molecular mechanisms and functioning of a variety of chemical and biomolecular systems. This includes the self-assembly and conformational stability of synthetic organic rosette nanotubes (RNTs), the aggregation of peptides and proteins related to neurodegeneration, the binding of ligands to proteins, and the solvation properties of biomolecules related to their functions. The replica RISM-KH-VM molecular solvation theory predicts and explains the structure, thermodynamics, and electrochemistry of electrolyte solutions sorbed in nanoporous carbon supercapacitor electrodes, and is part of recent research and development efforts. A new quasidynamics protocol couples multiple time step molecular dynamics (MTS-MD) stabilized with an optimized isokinetic Nosé–Hoover (OIN) thermostat driven by 3D-RISM-KH mean solvation forces at gigantic outer time steps of picoseconds, which are extrapolated forward at short inner time steps of femtoseconds with generalized solvation force extrapolation (GSFE). The OIN/3D-RISM-KH/GSFE quasidynamics is implemented in the Amber Molecular Dynamics package. It is validated on miniprotein 1L2Y and protein G in ambient aqueous solution, and shows the rate of sampling 150 times faster than in standard MD simulations on these biomolecules in explicit water. The self-consistent field version of Kohn–Sham DFT in 3D-RISM-KH mean solvation forces is implemented in the Amsterdam Density Functional (ADF) package. Its applications range from solvation thermochemistry, conformational equilibria, and photochemistry to activation barriers of different nanosystems in solutions and ionic liquids.

Protein pK a Prediction with Machine Learning.

Protein pKa prediction is essential for the investigation of the pH-associated relationship between protein structure and function. In this work, we introduce a deep learning-based protein pKa predictor DeepKa, which is trained and validated with the pKa values derived from continuous constant-pH molecular dynamics (CpHMD) simulations of 279 soluble proteins. Here, the CpHMD implemented in the Amber molecular dynamics package has been employed (HuangY.J. Chem. Inf. Model.2018, 58, 1372−138329949356). Notably, to avoid discontinuities at the boundary, grid charges are proposed to represent protein electrostatics. We show that the prediction accuracy by DeepKa is close to that by CpHMD benchmarking simulations, validating DeepKa as an efficient protein pKa predictor. In addition, the training and validation sets created in this study can be applied to the development of machine learning-based protein pKa predictors in the future. Finally, the grid charge representation is general and applicable to other topics, such as the protein–ligand binding affinity prediction.

Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery.

Predicting protein-ligand binding affinities and the associated thermodynamics of biomolecular recognition is a primary objective of structure-based drug design. Alchemical free energy simulations offer a highly accurate and computationally efficient route to achieving this goal. While the AMBER molecular dynamics package has successfully been used for alchemical free energy simulations in academic research groups for decades, widespread impact in industrial drug discovery settings has been minimal because of the previous limitations within the AMBER alchemical code, coupled with challenges in system setup and postprocessing workflows. Through a close academia-industry collaboration we have addressed many of the previous limitations with an aim to improve accuracy, efficiency, and robustness of alchemical binding free energy simulations in industrial drug discovery applications. Here, we highlight some of the recent advances in AMBER20 with a focus on alchemical binding free energy (BFE) calculations, which are less computationally intensive than alternative binding free energy methods where full binding/unbinding paths are explored. In addition to scientific and technical advances in AMBER20, we also describe the essential practical aspects associated with running relative alchemical BFE calculations, along with recommendations for best practices, highlighting the importance not only of the alchemical simulation code but also the auxiliary functionalities and expertise required to obtain accurate and reliable results. This work is intended to provide a contemporary overview of the scientific, technical, and practical issues associated with running relative BFE simulations in AMBER20, with a focus on real-world drug discovery applications.

4094 Structural Determinants of Immunogenicity for Peptide-Based Immunotherapy

OBJECTIVES/GOALS: Neoantigen vaccine immunotherapies have shown promise in clinical trials, but identifying which peptides to include in a vaccine remains a challenge. We aim to establish that molecular structural features can help predict which neoantigens to target to achieve tumor regression. METHODS/STUDY POPULATION: Proteins were prepared by recombinant expression in E. coli followed by in vitro refolding. Correctly folded proteins were purified by chromatography. Affinities of protein-protein interactions were measured by surface plasmon resonance (SPR) and thermal stabilities of proteins were determined by differential scanning fluorimetry. All experiments were performed at least in triplicate. Protein crystals were obtained by hanging drop vapor diffusion. The protein crystal structures were solved by molecular replacement and underwent several rounds of automated refinement. Molecular dynamics simulations were performed using the AMBER molecular dynamics package. RESULTS/ANTICIPATED RESULTS: A T cell receptor (TCR) expressed by tumor-infiltrating T cells exhibited a 20-fold stronger binding affinity to the neoantigen peptide compared to the self-peptide. X-ray crystal structures of the peptides with the major histocompatibility complex (MHC) protein demonstrated that a non-mutated residue in the peptide samples different positions with the mutation. The difference in conformations of the non-mutated residue was supported by molecular dynamics simulations. Crystal structures of the TCR engaging both peptide/MHCs suggested that the conformation favored by the mutant peptide was crucial for TCR binding. The TCR bound the neoantigen/MHC with faster binding kinetics. DISCUSSION/SIGNIFICANCE OF IMPACT: Our results suggest that the mutation impacts the conformation of another residue in the peptide, and this alteration allows for more favorable T cell receptor binding to the neoantigen. This highlights the potential of non-mutated residues in contributing to neoantigen recognition.

GPU-Accelerated Implementation of Continuous Constant pH Molecular Dynamics in Amber: pKa Predictions with Single-pH Simulations.

We present a GPU implementation of the continuous constant pH molecular dynamics (CpHMD) based on the most recent generalized Born implicit-solvent model in the pmemd engine of the Amber molecular dynamics package. To test the accuracy of the tool for rapid pKa predictions, a series of 2 ns single-pH simulations were performed for over 120 titratable residues in 10 benchmark proteins that were previously used to test the various continuous CpHMD methods. The calculated pKa's showed a root-mean-square deviation of 0.80 and correlation coefficient of 0.83 with respect to experiment. Also, 90% of the pKa's were converged with estimated errors below 0.1 pH units. Surprisingly, this level of accuracy is similar to our previous replica-exchange simulations with 2 ns per replica and an exchange attempt frequency of 2 ps-1 (Huang, Harris, and Shen J. Chem. Inf. Model. 2018 , 58 , 1372 - 1383 ). Interestingly, for the linked titration sites in two enzymes, although residue-specific protonation state sampling in the single-pH simulations was not converged within 2 ns, the protonation fraction of the linked residues appeared to be largely converged, and the experimental macroscopic pKa values were reproduced to within 1 pH unit. Comparison with replica-exchange simulations with different exchange attempt frequencies showed that the splitting between the two macroscopic pKa's is underestimated with frequent exchange attempts such as 2 ps-1, while single-pH simulations overestimate the splitting. The same trend is seen for the single-pH vs replica-exchange simulations of a hydrogen-bonded aspartyl dyad in a much larger protein. A 2 ns single-pH simulation of a 400-residue protein takes about 1 h on a single NVIDIA GeForce RTX 2080 graphics card, which is over 1000 times faster than a CpHMD run on a single CPU core of a high-performance computing cluster node. Thus, we envision that GPU-accelerated continuous CpHMD may be used in routine pKa predictions for a variety of applications, from assisting MD simulations with protonation state assignment to offering pH-dependent corrections of binding free energies and identifying reactive hot spots for covalent drug design.

Generalized Born Based Continuous Constant pH Molecular Dynamics in Amber: Implementation, Benchmarking and Analysis.

Solution pH plays an important role in structure and dynamics of biomolecular systems; however, pH effects cannot be accurately accounted for in conventional molecular dynamics simulations based on fixed protonation states. Continuous constant pH molecular dynamics (CpHMD) based on the λ-dynamics framework calculates protonation states on the fly during dynamical simulation at a specified pH condition. Here we report the CPU-based implementation of the CpHMD method based on the GBNeck2 generalized Born (GB) implicit-solvent model in the pmemd engine of the Amber molecular dynamics package. The performance of the method was tested using pH replica-exchange titration simulations of Asp, Glu and His side chains in 4 miniproteins and 7 enzymes with experimentally known p Ka's, some of which are significantly shifted from the model values. The added computational cost due to CpHMD titration ranges from 11 to 33% for the data set and scales roughly linearly as the ratio between the titrable sites and number of solute atoms. Comparison of the experimental and calculated p Ka's using 2 ns per replica sampling yielded a mean unsigned error of 0.70, a root-mean-squared error of 0.91, and a linear correlation coefficient of 0.79. Though this level of accuracy is similar to the GBSW-based CpHMD in CHARMM, in contrast to the latter, the current implementation was able to reproduce the experimental orders of the p Ka's of the coupled carboxylic dyads. We quantified the sampling errors, which revealed that prolonged simulation is needed to converge p Ka's of several titratable groups involved in salt-bridge-like interactions or deeply buried in the protein interior. Our benchmark data demonstrate that GBNeck2-CpHMD is an attractive tool for protein p Ka predictions.

Fast and flexible gpu accelerated binding free energy calculations within the amber molecular dynamics package

Alchemical free energy (AFE) calculations based on molecular dynamics (MD) simulations are key tools in both improving our understanding of a wide variety of biological processes and accelerating the design and optimization of therapeutics for numerous diseases. Computing power and theory have, however, long been insufficient to enable AFE calculations to be routinely applied in early stage drug discovery. One of the major difficulties in performing AFE calculations is the length of time required for calculations to converge to an ensemble average. CPU implementations of MD-based free energy algorithms can effectively only reach tens of nanoseconds per day for systems on the order of 50,000 atoms, even running on massively parallel supercomputers. Therefore, converged free energy calculations on large numbers of potential lead compounds are often untenable, preventing researchers from gaining crucial insight into molecular recognition, potential druggability and other crucial areas of interest. Graphics Processing Units (GPUs) can help address this. We present here a seamless GPU implementation, within the PMEMD module of the AMBER molecular dynamics package, of thermodynamic integration (TI) capable of reaching speeds of >140 ns/day for a 44,907-atom system, with accuracy equivalent to the existing CPU implementation in AMBER. The implementation described here is currently part of the AMBER 18 beta code and will be an integral part of the upcoming version 18 release of AMBER. © 2018 Wiley Periodicals, Inc.

Protocols for Molecular Dynamics Simulations of RNA Nanostructures.

Molecular dynamics (MD) simulations have been used as one of the main research tools to study a wide range of biological systems and bridge the gap between X-ray crystallography or NMR structures and biological mechanism. In the field of RNA nanostructures, MD simulations have been used to fix steric clashes in computationally designed RNA nanostructures, characterize the dynamics, and investigate the interaction between RNA and other biomolecules such as delivery agents and membranes.In this chapter we present examples of computational protocols for molecular dynamics simulations in explicit and implicit solvent using the Amber Molecular Dynamics Package. We also show examples of post-simulation analysis steps and briefly mention selected tools beyond the Amber package. Limitations of the methods, tools, and protocols are also discussed. Most of the examples are illustrated for a small RNA duplex (helix), but the protocols are applicable to any nucleic acid structure, subject only to the computational speed and memory limitations of the hardware available to the user.

Exploiting molecular dynamics in Nested Sampling simulations of small peptides

Nested Sampling (NS) is a parameter space sampling algorithm which can be used for sampling the equilibrium thermodynamics of atomistic systems. NS has previously been used to explore the potential energy surface of a coarse-grained protein model and has significantly outperformed parallel tempering when calculating heat capacity curves of Lennard-Jones clusters. The original NS algorithm uses Monte Carlo (MC) moves; however, a variant, Galilean NS, has recently been introduced which allows NS to be incorporated into a molecular dynamics framework, so NS can be used for systems which lack efficient prescribed MC moves. In this work we demonstrate the applicability of Galilean NS to atomistic systems. We present an implementation of Galilean NS using the Amber molecular dynamics package and demonstrate its viability by sampling alanine dipeptide, both in vacuo and implicit solvent. Unlike previous studies of this system, we present the heat capacity curves of alanine dipeptide, whose calculation provides a stringent test for sampling algorithms. We also compare our results with those calculated using replica exchange molecular dynamics (REMD) and find good agreement. We show the computational effort required for accurate heat capacity estimation for small peptides. We also calculate the alanine dipeptide Ramachandran free energy surface for a range of temperatures and use it to compare the results using the latest Amber force field with previous theoretical and experimental results.

Evaluating thermodynamic integration performance of the new amber molecular dynamics package and assess potential halogen bonds of enoyl-ACP reductase (FabI) benzimidazole inhibitors

Thermodynamic integration (TI) can provide accurate binding free energy insights in a lead optimization program, but its high computational expense has limited its usage. In the effort of developing an efficient and accurate TI protocol for FabI inhibitors lead optimization program, we carefully compared TI with different Amber molecular dynamics (MD) engines (sander and pmemd), MD simulation lengths, the number of intermediate states and transformation steps, and the Lennard-Jones and Coulomb Softcore potentials parameters in the one-step TI, using eleven benzimidazole inhibitors in complex with Francisella tularensis enoyl acyl reductase (FtFabI). To our knowledge, this is the first study to extensively test the new AMBER MD engine, pmemd, on TI and compare the parameters of the Softcore potentials in the one-step TI in a protein-ligand binding system. The best performing model, the one-step pmemd TI, using 6 intermediate states and 1 ns MD simulations, provides better agreement with experimental results (RMSD = 0.52 kcal/mol) than the best performing implicit solvent method, QM/MM-GBSA from our previous study (RMSD = 3.00 kcal/mol), while maintaining similar efficiency. Briefly, we show the optimized TI protocol to be highly accurate and affordable for the FtFabI system. This approach can be implemented in a larger scale benzimidazole scaffold lead optimization against FtFabI. Lastly, the TI results here also provide structure-activity relationship insights, and suggest the parahalogen in benzimidazole compounds might form a weak halogen bond with FabI, which is a well-known halogen bond favoring enzyme.

A Molecular Dynamics Study of Interaction between Kinesin Light Chain and Amyloid Precursor Protein

The misregulation of the cargo binding domain of kinesin‐1 light chain (KLC) has been linked to various neurodegenerative diseases, most notably Alzheimer's, as KLC may play a role in amyloid precursor protein (APP) transport and subsequent degradation. Recent literature suggests a third protein, C‐Jun‐amino‐terminal kinase‐interacting protein 1 (JIP1), acts as a mediator protein for binding of KLC to APP. Computational molecular dynamics (MD) was conducted on this ternary complex, as well as the JIP1:KLC complex, using the Amber MD package to understand the atomistic interactions between the structures. Schrodinger's Glide was used to dock the structure of JIP1 (extracted from PDB 2H96) to KLC to form the binary complex, while the online Patchdock server was used to dock APP to the JIP1:KLC complex to form the ternary state. Within simulations of various JIP1:KLC complexes, the JIP1 protein would orient itself in nearly the same location relative to KLC regardless of the initial starting conformation used for MD. The simulations were also used to help elucidate a low energy conformation of the bound complex, since this conformation has not been resolved experimentally. Preliminary binding free energy calculations suggest that the ternary complex is more stable than the binary complex (APP:KLC).

Computational Analysis of the Interaction between Kinesin Light Chain and Amyloid Precursor Protein

Disruption of binding between the cargo binding domain of kinesin‐1 light chain (KLC) and amyloid precursor protein (APP) has been implicated in a variety of diseases, including familial Alzheimer's. KLC serves to transport cargo, such as APP, across a microtubule; however, when this function is disrupted, APP releases and disrupts neuronal function. In this project, computational molecular dynamics (MD) was conducted on apo KLC, apo APP, and docked KLC‐APP structures using the Amber MD package to understand the dynamics and interactions between the structures in atomistic detail. KLC remained stable in the apo state and helped to stabilize APP when docked. Preliminary per‐residue free energy decomposition was utilized to determine the key residues involved in KLC‐APP binding. MD was also conducted on an extended tetratricopeptide repeat (TPR) structure to observe fast‐folding dynamics and produce the correct alpha‐helical formation. TPR units, in tandem, form KLC's cargo‐binding domain. Within 700 ns, TPR was able to fold to mostly resemble the consensus crystal structure (CTPR) to within an RMSD of 2.13 Å.

Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald.

We present an implementation of explicit solvent all atom classical molecular dynamics (MD) within the AMBER program package that runs entirely on CUDA-enabled GPUs. First released publicly in April 2010 as part of version 11 of the AMBER MD package and further improved and optimized over the last two years, this implementation supports the three most widely used statistical mechanical ensembles (NVE, NVT, and NPT), uses particle mesh Ewald (PME) for the long-range electrostatics, and runs entirely on CUDA-enabled NVIDIA graphics processing units (GPUs), providing results that are statistically indistinguishable from the traditional CPU version of the software and with performance that exceeds that achievable by the CPU version of AMBER software running on all conventional CPU-based clusters and supercomputers. We briefly discuss three different precision models developed specifically for this work (SPDP, SPFP, and DPDP) and highlight the technical details of the approach as it extends beyond previously reported work [Götz et al., J. Chem. Theory Comput. 2012, DOI: 10.1021/ct200909j; Le Grand et al., Comp. Phys. Comm. 2013, DOI: 10.1016/j.cpc.2012.09.022].We highlight the substantial improvements in performance that are seen over traditional CPU-only machines and provide validation of our implementation and precision models. We also provide evidence supporting our decision to deprecate the previously described fully single precision (SPSP) model from the latest release of the AMBER software package.

Improving the Efficiency of Free Energy Calculations in the Amber Molecular Dynamics Package

Alchemical transformations are widely used methods to calculate free energies. Amber has traditionally included support for alchemical transformations as part of the sander molecular dynamics (MD) engine. Here we describe the implementation of a more efficient approach to alchemical transformations in the Amber MD package. Specifically we have implemented this new approach within the more computational efficient and scalable pmemd MD engine that is included with the Amber MD package. The majority of the gain in efficiency comes from the improved design of the calculation, which includes better parallel scaling and reduction in the calculation of redundant terms. This new implementation is able to reproduce results from equivalent simulations run with the existing functionality, but at 2.5 times greater computational efficiency. This new implementation is also able to run softcore simulations at the λ end states making direct calculation of free energies more accurate, compared to the extrapolation required in the existing implementation. The updated alchemical transformation functionality will be included in the next major release of Amber (scheduled for release in Q1 2014) and will be available at http://ambermd.org, under the Amber license.

Electronic Absorption Spectra from MM and ab initio QM/MM Molecular Dynamics: Environmental Effects on the Absorption Spectrum of Photoactive Yellow Protein.

We describe a new interface of the GPU parallelized TeraChem electronic structure package and the Amber molecular dynamics package for quantum mechanical (QM) and mixed QM and molecular mechanical (MM) molecular dynamics simulations. This QM/MM interface is used for computation of the absorption spectra of the photoactive yellow protein (PYP) chromophore in vacuum, aqueous solution, and protein environments. The computed excitation energies of PYP require a very large QM region (hundreds of atoms) covalently bonded to the chromophore in order to achieve agreement with calculations that treat the entire protein quantum mechanically. We also show that 40 or more surrounding water molecules must be included in the QM region in order to obtain converged excitation energies of the solvated PYP chromophore. These results indicate that large QM regions (with hundreds of atoms) are a necessity in QM/MM calculations.

Multiscale methods for nanochemistry and biophysics in solution

Statistical–mechanical, 3D-RISM-KH molecular theory of solvation is an integral equation theory of molecular liquids capable of describing, from the first principles, solvation structure and thermodynamics in complex nanoscale systems. The 3D-RISM-KH theory accurately yields the solvation structure for biomolecular systems as large as chaperonins and ion channels. It predicts 3D maps of ligand binding affinity at once without using empirical scoring functions, and has significant potential for computer-aided drug design. Recently, the 3D-RISM-KH theory was coupled with MD simulations by contracting solvent degrees of freedom to treat solvated biomolecules, and the MD/3D-RISM-KH multiscale method was implemented in the Amber molecular dynamics package [T. Luchko et al., 2010, [ 29]]. The MD/3D-RISM-KH method allows one to study biomolecular processes on extremely long timescales, as the statistics of rare solvent and ligand events is accounted for analytically. We further replaced the MM/GB(PB)SA post-processing which evaluates the free energy by empirical treatment of non-polar solvation effects with the MM/3D-RISM-KH method which employs 3D-RISM-KH to evaluate the solvation thermodynamics. In this paper, we briefly review the 3D-RISM-KH theory as well as the MD/3D-RISM-KH multiscale method implemented in the Amber package. We illustrate the MM/3D-RISM-KH method on binding of thiamine (vitamin of group B) to extracytoplasmic thiamine-binding lipoprotein MG289, and on the microscopic solvation properties of the bacterial Gloeobacter violaceus pentameric ligand-gated ion channel (GLIC) homologue.

Implementation of the SCC-DFTB method for hybrid QM/MM simulations within the amber molecular dynamics package.

Self-consistent charge density functional tight-binding (SCC-DFTB) is a semiempirical method based on density functional theory and has in many cases been shown to provide relative energies and geometries comparable in accuracy to full DFT or ab initio MP2 calculations using large basis sets. This article shows an implementation of the SCC-DFTB method as part of the new QM/MM support in the AMBER 9 molecular dynamics program suite. Details of the implementation and examples of applications are shown.

Swelling behavior of the cellulose Iβ crystal models by molecular dynamics

The various crystal models of cellulose Iβ, each differing in crystal size, have been studied by computer simulation using the amber molecular-dynamics package and the GLYCAM parameters. The four types of crystal model were constructed by a combination of two base-plane sizes, consisting of either 24 or 48 chains and two chain lengths having either 10 or 20 residues. The base planes of the crystal models were composed by the edges of the [1, 1, 0], [1, −1, 0], and [1, 0, 0] crystal planes, where the [1, 1, 0] plane was assigned to the longest edge. The crystal models were soaked in water boxes to investigate their swelling behavior. Unexpectedly, the crystal models twisted quickly to form a slightly right-handed shape during the initial ∼50 ps and that, in a steady, swollen state, the twisted forms remained for the rest of the simulation time. In spite of such overall deformation, the inner part of the swollen model fairly reproduced the important structural features of the original crystal structure, such as the rotational positions of the substituent groups and the hydrogen-bonding scheme. On heating the crystal model up to 550 K, the twisted shape was conserved in most of the temperature range, while the initial conformations of the substituent groups deviated above ∼430 K, followed by appreciable disordering in chain sheets at higher temperatures. It is suggested that some internal tensions are involved within a chain sheet of the initial structure. In the course of swelling, some of these tensions were released to introduce a twisted shape in the crystal models.