Standard Molecular Dynamics Research Articles

BLaDE: A Basic Lambda Dynamics Engine for GPU-Accelerated Molecular Dynamics Free Energy Calculations.

There is an accelerating interest in practical applications of alchemical free energy methods to problems in protein design, constant pH simulations, and especially computer-aided drug design. In the present paper, we describe a basic lambda dynamics engine (BLaDE) that enables alchemical free energy simulations, including multisite λ dynamics (MSλD) simulations, on graphical processor units (GPUs). We find that BLaDE is 5 to 8 times faster than the current GPU implementation of MSλD-based free energy calculations in CHARMM. We also demonstrate that BLaDE running standard molecular dynamics attains a performance competitive with and sometimes exceeding that of the highly optimized OpenMM GPU code. BLaDE is available as a standalone program and through an API in CHARMM.

Hybrid polarizable simulations of a conventional hydrophobic polyelectrolyte. Toward a theoretical tool for green science innovation.

Hybrid modeling approaches based on all-atom force fields to handle a solute and coarse-grained models to account for the solvent are promising numerical tools that can be used to understand the properties of large and multi-components solutions and thus to speed up the development of new industrial products that obey the standard of green and sustainable chemistry. Here, we discuss the ability of a full polarizable hybrid approach coupled to a standard molecular dynamics scheme to model the behavior in the aqueous phase and at infinite dilution conditions of a standard hydrophobic polyelectrolyte polymer whose charge is neutralized by explicit counterions. Beyond the standard picture of a polyelectrolyte behavior governed by an interplay between opposite intra-polyelectrolyte and inter-polyelectrolyte/counterion Coulombic effects, our simulations show the key role played by both intra-solute polarization effects and long range solute/solvent electrostatics to stabilize compact globular conformations of that polyelectrolyte. Our full polarizable hybrid modeling approach is thus a new theoretical tool well suited to be used in digital strategies for accelerating innovation for green science, for instance.

Simulation of protein pulling dynamics on second time scale with boxed molecular dynamics.

We demonstrate how recently developed Boxed Molecular Dynamics (BXD) and kinetics [D. V. Shalashilin et al., J. Chem. Phys. 137, 165102 (2012)] can provide a kinetic description of protein pulling experiments, allowing for a connection to be made between experiment and the atomistic protein structure. BXD theory applied to atomic force microscopy unfolding is similar in spirit to the kinetic two-state model [A. Noy and R. W. Friddle, Methods 60, 142 (2013)] but with some differences. First, BXD uses a large number of boxes, and therefore, it is not a two-state model. Second, BXD rate coefficients are obtained from atomistic molecular dynamics simulations. BXD can describe the dependence of the pulling force on pulling speed. Similar to Shalashilin et al. [J. Chem. Phys. 137, 165102 (2012)], we show that BXD is able to model the experiment at a very long time scale up to seconds, which is way out of reach for standard molecular dynamics.

Physical Origin of Distinct Mechanical Properties of Polymer Tethered Graphene Nanosheets Reinforced Polymer Nanocomposites Revealed by Nonequilibrium Molecular Dynamics Simulations

AbstractA non‐equilibrium deformation simulation based on standard molecular dynamics is employed to reveal the physical origin of distinct mechanical properties of polymer nanocomposites (PNCs) reinforced by polymer tethered graphene nanosheets (P/G‐T). The effect of tethered polymer length (L) and interaction strength on the mechanical properties, including tensile stress, modulus, and yield strength, are examined. The simulation results show that the P/G‐T systems exhibit improved mechanical properties as L and interaction strength increases. The strengthening in attractive interaction between matrix polymer and tethered polymer has a better effect on mechanical properties than that between P and graphene nanosheets. It is found that the bond orientation and nonbonding potential is of crucial importance in determining the mechanical properties of the P/G‐T nanocomposite systems. The calculations of radial distribution functions and mean‐squared displacement are further performed to reveal the physical origin of enhanced mechanical properties, suggesting that the stronger interfacial interactions will induce the closer packing distance with smaller free volume, higher polymer chain entanglement, and higher restriction of polymer chain movement, which synergistically contribute to the improvement in the mechanical properties. The results may provide useful guidance for promoting the development and practical application of the advanced PNCs with excellent performance.

Exploring the conformational space of a receptor for drug design: An ERα case study

Protein flexibility is challenging for both experimentalists and modellers, especially in the field of drug design. Estrogen Receptor alpha (ERα) is an extensively studied Nuclear Receptor (NR) and a well-known therapeutic target with an important role in development and physiology. It is also a frequent off-target in standard toxicity tests for endocrine disruption. Here, we aim to evaluate the degree to which the conformational space and macromolecular flexibility of this well-characterized drug target can be described. Our approach exploits hundreds of crystallographic structures by means of molecular dynamics simulations and of virtual screening. The analysis of hundreds of crystal structures confirms the presence of two main conformational states, known as ’agonist’ and ’antagonist’, that mainly differ by the orientation of the C-terminal helix H12 which serves to close the binding pocket. ERα also shows some loop flexibility, as well as variable side-chain orientations in its active site. We scrutinized the extent to which standard molecular dynamics simulations or crystallographic refinement as ensemble recapitulate most of the variability features seen by the array of available crystal structures. In parallel, we investigated on the kind and extent of flexibility that are required to achieve convincing docking for all high-affinity ERα ligands present in BindingDB. Using either only one conformation with a few side-chains set flexible, or static structure ensembles in parallel during docking led to good quality and similar pose predictions. These results suggest that the several hundreds of crystal structures already known can properly describe the whole conformational universe of ERα′s ligand binding domain. This opens the road for better drug design and affinity computation.

(Invited) Computational Simulations of Selective Interactions between ssDNA and SWCNTs

Structure-specific interactions between single-stranded DNA oligos and SWCNTs form the basis of important selective sorting methods yet remain poorly understood. Among these interactions is the differential affinity of (ATT)4 for one of the two (7,5) enantiomers. To understand the nature of the affinity difference, we have performed molecular dynamics (MD) simulations followed by free energy (replica exchange) calculations over the temperature range of 290 to 727 K. Comparisons between the conformational free energy landscapes of (ATT)4 structure on right-handed and left-handed (7,5) SWCNTs (respectively known as M and P in stereochemical convention) suggest that (ATT)4 has more contact area with M than with the P enantiomer. However, the corresponding end-to-end distances of (ATT)4 fall nearly in the same region on the landscape. Our standard MD simulations for (ATTT)3 interactions with SWCNTs also suggest that this DNA sequence covers (8,4) SWCNT more effectively than (7,6) and (8,3) SWCNTs, which is consistent with our recent dye-quenching experiments.Other experiments have revealed that when guanines in ssDNA oligos are covalently functionalized to SWCNT sidewalls, the extent of functionalization can depend on the position of the guanine within the oligo. Central positions within long oligos (e.g. T15GT15) lead to greater functionalization than end guanine positions (e.g. T30G). We have performed both standard MD and free energy calculations to investigate the source of this effect. For efficient guanine functionalization of SWCNTs, guanine needs to remain adsorbed on the SWCNT surface. Our computational observations suggest that the central guanine within T15GT15 oligo remains adsorbed on the (6,5) SWCNT surface even when the DNA interacts with other neighboring DNA strands. Under these conditions, the end guanine in T30G is destabilized and tends to desorb from (6,5) SWCNT surface due to inter-strand interactions. Our free energy calculations also suggest that oligos with a central guanine assume more elongated DNA conformations at the surface of a (6,5) SWCNT than do oligos with end guanines. We suggest that DNA elongation is associated with more stable adsorption on the SWCNT surface, supporting efficient covalent functionalization of central guanines. Figure 1

Particle-particle particle-mesh algorithm for electrolytes between charged dielectric interfaces.

Ionic distributions near charged interfaces control processes from colloidal aggregation to solvent flow in nanodevices. Such interfaces are often characterized by a jump in the permittivity, which gives rise to the surface polarization charge. This induced charge may significantly affect the ionic distributions so that efficient methods for modeling spatially varying dielectrics are needed. We formulate a method with O(N⁡log⁡N) scaling for electrolytes between charged planar interfaces with asymmetric dielectric contrasts. Our approach, which builds on earlier work, is based on combining image charges with the particle-particle particle-mesh algorithm and representing uniform surface charges via an electric field. This enables simulations of complex dielectric interactions that outperform most alternative methods in speed and accuracy. To make the method practically useful, we provide guidelines-based upon careful tests-for choosing optimal simulation parameters. Explicit expressions for the electrostatic forces are given to facilitate the implementation of our algorithm in standard molecular dynamics packages.

Spike protein fusion loop controls SARS-CoV-2 fusogenicity and infectivity

The high SARS-CoV-2 reproductive number driving the COVID-19 pandemic has been a mystery. Our recent in vitro, and in vivo coronaviral pathogenesis studies involving Mouse Hepatitis Virus (MHV-A59) suggest a crucial role for a small host membrane-virus contact initiator region of the Spike protein, called the fusion peptide that enhances the virus fusogenicity and infectivity. Here I study the Spike from five human β-coronaviruses (HCoV) including the SARS-CoV-2, and MHV-A59 for comparison. The structural and dynamics analyses of the Spike show that its fusion loop spatially organizes three fusion peptides contiguous to each other to synergistically trigger the virus-host membrane fusion process. I propose a Contact Initiation Model based on the architecture of the Spike quaternary structure that explains the obligatory participation of the fusion loop in the initiation of the host membrane contact for the virus fusion process. Among all the HCoV Spikes in this study, SARS-CoV-2 has the most hydrophobic surface and the extent of hydrophobicity correlates with the reproductive number and infectivity of the other HCoV. Comparison between results from standard and replica exchange molecular dynamics reveal the unique physicochemical properties of the SARS-CoV-2 fusion peptides, accrued in part from the presence of consecutive prolines that impart backbone rigidity which aids the virus fusogenicity. The priming of the Spike by its cleavage and subsequent fusogenic conformational transition steered by the fusion loop may be critical for the SARS-CoV-2 spread. The importance of the fusion loop makes it an apt target for anti-virals and vaccine candidates.

Full structural ensembles of intrinsically disordered proteins from unbiased molecular dynamics simulations

Molecular dynamics (MD) simulation is widely used to complement ensemble-averaged experiments of intrinsically disordered proteins (IDPs). However, MD often suffers from limitations of inaccuracy. Here, we show that enhancing the sampling using Hamiltonian replica-exchange MD (HREMD) led to unbiased and accurate ensembles, reproducing small-angle scattering and NMR chemical shift experiments, for three IDPs of varying sequence properties using two recently optimized force fields, indicating the general applicability of HREMD for IDPs. We further demonstrate that, unlike HREMD, standard MD can reproduce experimental NMR chemical shifts, but not small-angle scattering data, suggesting chemical shifts are insufficient for testing the validity of IDP ensembles. Surprisingly, we reveal that despite differences in their sequence, the inter-chain statistics of all three IDPs are similar for short contour lengths (< 10 residues). The results suggest that the major hurdle of generating an accurate unbiased ensemble for IDPs has now been largely overcome.

Energy-dependent protein folding: modeling how a protein folding machine may work

Background: Proteins fold robustly and reproducibly in vivo, but many cannot fold in vitro in isolation from cellular components. Despite the remarkable progress that has been achieved by the artificial intelligence approaches in predicting the protein native conformations, the pathways that lead to such conformations, either in vitro or in vivo, remain largely unknown. The slow progress in recapitulating protein folding pathways in silico may be an indication of the fundamental deficiencies in our understanding of folding as it occurs in nature. Here we consider the possibility that protein folding in living cells may not be driven solely by the decrease in Gibbs free energy and propose that protein folding in vivo should be modeled as an active energy-dependent process. The mechanism of action of such a protein folding machine might include direct manipulation of the peptide backbone. Methods: To show the feasibility of a protein folding machine, we conducted molecular dynamics simulations that were augmented by the application of mechanical force to rotate the C-terminal amino acid while simultaneously limiting the N-terminal amino acid movements. Results: Remarkably, the addition of this simple manipulation of peptide backbones to the standard molecular dynamics simulation indeed facilitated the formation of native structures in five diverse alpha-helical peptides. Steric clashes that arise in the peptides due to the forced directional rotation resulted in the behavior of the peptide backbone no longer resembling a freely jointed chain. Conclusions: These simulations show the feasibility of a protein folding machine operating under the conditions when the movements of the polypeptide backbone are restricted by applying external forces and constraints. Further investigation is needed to see whether such an effect may play a role during co-translational protein folding in vivo and how it can be utilized to facilitate folding of proteins in artificial environments.

An algorithm for estimating kinetic parameters of atomistic rare events using finite-time temperature programmed molecular dynamics trajectories

A new variation of the temperature programmed molecular dynamics (TPMD) algorithm, called finite-time TPMD (FT-TPMD), is presented. In this variation, a collection of independent molecular dynamics (MD) trajectories is generated, such that each trajectory is of a fixed duration of time and a transition from a given state of the system may/may not occur in a trajectory. This eliminates the requirement of a successful transition in each trajectory and on-the-fly detection of transitions needed in previous TPMD implementations. As a result, FT-TPMD can be easily implemented with standard MD packages, such as LAMMPS, for performing TPMD calculations. Compared to previous TPMD algorithm, FT-TPMD is shown to be computationally as efficient in some cases and better in others. In addition, several refinements are discussed that enable systematic detection of rare transition events, classification of states, and calculation of rate constants. A set of rules is devised to identify one or more atoms involved in a transition and disqualify recrossing events.

Insights into the Polymorphic Structures and Enantiotropic Layer-Slip Transition in Paracetamol Form III from Enhanced Molecular Dynamics

Reversible temperature-mediated solid phase changes, otherwise known as enantiotropic transformations, occur in many molecular crystals. These transformations take place as a result of the free-energy stabilization through entropic contributions at finite temperatures and can often have significant implications for the properties of crystalline solids. As such, understanding and predicting these transformations is of great importance. In this study, we utilize molecular simulations to elucidate the mechanism behind the enantiotropic layer-slip phase transformation between the orthorhombic and monoclinic versions of paracetamol form III (form III-m and III-o). Using standard molecular dynamics (MD) in addition to crystal adiabatic free-energy dynamics, an MD-based enhanced sampling approach for crystalline systems, we demonstrate that the transformation from the monoclinic form III-m to the orthorhombic form III-o is driven by localized and dynamic disorder within the structure rather than a perfect crystal-to-crystal transition. These results suggest that the orthorhombic form III-o structure does not exist as a perfect orthorhombic crystal with fully aligned layers, but rather, as an entropy-stabilized collective ensemble average of various misaligned layer-slipped structures. Overall, these simulation approaches, which explicitly treat dynamic structural disorder, allowed us to map out the free-energy landscape of this enantiotropic transformation as a function of temperature and extract critical insights into the underlying mechanism of the transformation.

Energy-dependent protein folding: modeling how a protein folding machine may work.

Background: Proteins fold robustly and reproducibly in vivo, but many cannot fold in vitro in isolation from cellular components. Despite the remarkable progress that has been achieved by the artificial intelligence approaches in predicting the protein native conformations, the pathways that lead to such conformations, either in vitro or in vivo, remain largely unknown. The slow progress in recapitulating protein folding pathways in silico may be an indication of the fundamental deficiencies in our understanding of folding as it occurs in nature. Here we consider the possibility that protein folding in living cells may not be driven solely by the decrease in Gibbs free energy and propose that protein folding in vivo should be modeled as an active energy-dependent process. The mechanism of action of such a protein folding machine might include direct manipulation of the peptide backbone. Methods: To show the feasibility of a protein folding machine, we conducted molecular dynamics simulations that were augmented by the application of mechanical force to rotate the C-terminal amino acid while simultaneously limiting the N-terminal amino acid movements. Results: Remarkably, the addition of this simple manipulation of peptide backbones to the standard molecular dynamics simulation indeed facilitated the formation of native structures in five diverse alpha-helical peptides. Steric clashes that arise in the peptides due to the forced directional rotation resulted in the behavior of the peptide backbone no longer resembling a freely jointed chain. Conclusions: These simulations show the feasibility of a protein folding machine operating under the conditions when the movements of the polypeptide backbone are restricted by applying external forces and constraints. Further investigation is needed to see whether such an effect may play a role during co-translational protein folding in vivo and how it can be utilized to facilitate folding of proteins in artificial environments.

Extended diffusion theory: Recovering dynamics from biased/accelerated molecular simulations.

Dynamical properties are of great importance in determining the behavior of synthetic and natural molecules, but capturing them by computational methods is a nontrivial task. Very often the time scales of the relevant phenomena are far beyond the typical time windows accessible by classical Molecular Dynamics (MD) simulations, currently limited to the order of microseconds on standard laboratory workstations. On the other hand, biased and accelerated simulations allow for fast and thorough exploration of the molecular conformational space, but they lose the dynamic information. The problem of recovering dynamics from biased/accelerated simulations is a very active field of research, but no totally robust/reliable solutions have been given yet. In this paper it is shown how the Smoluchowski equation, in the framework of Diffusion Theory (DT), can be used to bridge this gap, and dynamical properties, in the form of time correlation functions (TCFs), can be extracted also from such kind of simulations. DT is first extended (EDT) to express the mobility tensors entering the Smoluchowski operator in terms of a recently introduced unified and regularized Rotne-Prager-Yamakawa approximation, [P. J. Zuk, E. Wajnryb, K. A. Mizerski, P. Szymczak, J. Fluid. Mech. 2014, 741, R5, 1-13] also involving mixed rotation-translation contributions, and rotation-rotation terms beside the classical translation-translation ones, so far used in DT. Then, the method is applied to recover the dynamics of a nontrivial example of a peptide in explicit water from the first 200 ns of a Replica Exchange Molecular Dynamics simulation, which is a popular computational method that destroys the long time dynamics. EDT dynamics were found to favorably compare against those coming from a standard MD simulation of the same system, requiring a time window of 30 μs to converge. This result shows that EDT is a tool of practical value to recover the long time dynamics of systems in diffusive regimes from biased/accelerated simulations, to be exploited in those cases when direct evaluation by standard MD is unfeasible.

Accounting for the Central Role of Interfacial Water in Protein-Ligand Binding Free Energy Calculations.

Rigorous binding free energy methods in drug discovery are growing in popularity because of a combination of methodological advances, improvements in computer hardware, and workflow automation. These calculations typically use molecular dynamics (MD) to sample from the Boltzmann distribution of conformational states. However, when part or all of the binding sites is inaccessible to the bulk solvent, the time needed for water molecules to equilibrate between bulk solvent and the binding site can be well beyond what is practical with standard MD. This sampling limitation is problematic in relative binding free energy calculations, which compute the reversible work of converting ligand 1 to ligand 2 within the binding site. Thus, if ligand 1 is smaller and/or more polar than ligand 2, the perturbation may allow additional water molecules to occupy a region of the binding site. However, this change in hydration may not be captured by standard MD simulations and may therefore lead to errors in the computed free energy. We recently developed a hybrid Monte Carlo/MD (MC/MD) method, which speeds up the equilibration of water between bulk solvent and buried cavities, while sampling from the intended distribution of states. Here, we report on the use of this approach in the context of alchemical binding free energy calculations. We find that using MC/MD markedly improves the accuracy of the calculations and also reduces hysteresis between the forward and reverse perturbations, relative to matched calculations using only MD with or without the crystallographic water molecules. The present method is available for use in AMBER simulation software.

Extensible and Scalable Adaptive Sampling on Supercomputers.

The accurate sampling of protein dynamics is an ongoing challenge despite the utilization of high-performance computer (HPC) systems. Utilizing only "brute force" molecular dynamics (MD) simulations requires an unacceptably long time to solution. Adaptive sampling methods allow a more effective sampling of protein dynamics than standard MD simulations. Depending on the restarting strategy, the speed up can be more than 1 order of magnitude. One challenge limiting the utilization of adaptive sampling by domain experts is the relatively high complexity of efficiently running adaptive sampling on HPC systems. We discuss how the ExTASY framework can set up new adaptive sampling strategies and reliably execute resulting workflows at scale on HPC platforms. Here, the folding dynamics of four proteins are predicted with no a priori information.

Insight into the Dual inhibitory Mechanism of verbascoside targeting serine/threonine phosphatase Stp1 against Staphylococcus aureus

The eukaryotic-like serine/threonine phosphatase (Stp1) is an enzyme-dependent protein phosphatase involved in regulating various virulence factors of Staphylococcus aureus. Owing to its role in S. aureus infections, Stp1 has become a potential target for antibiotic development. Unfortunately, there are very few reports describing Stp1 inhibitors. Using virtual screening, we have identified a potent and effective Stp1 inhibitor, verbascoside (VBS). Interestingly, the kinetics of the enzymatic reaction revealed that this natural inhibitor acts via both competitive and allosteric mechanisms. To explore the mechanism of interaction between VBS and Stp1, standard molecular dynamics (MD) simulations were performed for the Stp1-VBS complex. Consistent with the experimental results, competitive and allosteric binding sites for VBS were identified in Stp1. Met39, Gly41, His42, Arg161, and Asn162 residues were involved in the competitive binding of VBS, while Arg122, Ser136, Asp137, Asn142, and Val145 residues were associated with the allosteric binding of VBS. The contributions of these residues were confirmed by amino acid site-directed mutagenesis and fluorescence quenching experiments. This work demonstrates that VBS is a potent anti-virulence compound against S. aureus infection, laying the foundation for the further development of novel anti-virulence agents.

Casimir force and its effects on pull-in instability modelled using molecular dynamics simulations.

We present a new methodology to incorporate the Casimir forces within the molecular dynamics (MD) framework. At atomistic scales, the potential energy between two particles arising due to the Casimir effect can be represented as U(r ij ) = C/r 7. Incorporating the Casimir effect in MD simulations requires the knowledge of C, a problem hitherto unsolved. We overcome this by equating the total potential energy contributions due to each atomistic pair with the potential energy of continuum scale interacting bodies having similar geometries. After having identified the functional form of C, standard MD simulations are augmented with the potential energy contribution due to pairwise Casimir interactions. The developed framework is used to study effects of the Casimir force on the pull-in instability of rectangular and hollow cylindrical shaped deformable electrodes separated by a small distance from a fixed substrate electrode. Our MD results for pull-instability qualitatively agree with the previously reported analytical results but are quantitatively different. The effect of using longer-ranged Casimir forces in a constant temperature environment on the pull-in behaviour has also been studied.

A combination of machine learning and infrequent metadynamics to efficiently predict kinetic rates, transition states, and molecular determinants of drug dissociation from G protein-coupled receptors.

Determining the drug-target residence time (RT) is of major interest in drug discovery given that this kinetic parameter often represents a better indicator of in vivo drug efficacy than binding affinity. However, obtaining drug-target unbinding rates poses significant challenges, both computationally and experimentally. This is particularly palpable for complex systems like G Protein-Coupled Receptors (GPCRs) whose ligand unbinding typically requires very long timescales oftentimes inaccessible by standard molecular dynamics simulations. Enhanced sampling methods offer a useful alternative, and their efficiency can be further improved by using machine learning tools to identify optimal reaction coordinates. Here, we test the combination of two machine learning techniques, automatic mutual information noise omission and reweighted autoencoded variational Bayes for enhanced sampling, with infrequent metadynamics to efficiently study the unbinding kinetics of two classical drugs with different RTs in a prototypic GPCR, the μ-opioid receptor. Dissociation rates derived from these computations are within one order of magnitude from experimental values. We also use the simulation data to uncover the dissociation mechanisms of these drugs, shedding light on the structures of rate-limiting transition states, which, alongside metastable poses, are difficult to obtain experimentally but important to visualize when designing drugs with a desired kinetic profile.

Parallel algorithms for hyperdynamics and local hyperdynamics.

Hyperdynamics (HD) is a method for accelerating the timescale of standard molecular dynamics (MD). It can be used for simulations of systems with an energy potential landscape that is a collection of basins, separated by barriers, where transitions between basins are infrequent. HD enables the system to escape from a basin more quickly while enabling a statistically accurate renormalization of the simulation time, thus effectively boosting the timescale of the simulation. In the work of Kim et al. [J. Chem. Phys. 139, 144110 (2013)], a local version of HD was formulated, which exploits the intrinsic locality characteristic typical of most systems to mitigate the poor scaling properties of standard HD as the system size is increased. Here, we discuss how both HD and local HD can be formulated to run efficiently in parallel. We have implemented these ideas in the LAMMPS MD code, which means HD can be used with any interatomic potential LAMMPS supports. Together, these parallel methods allow simulations of any size to achieve the time acceleration offered by HD (which can be orders of magnitude), at a cost of 2-4× that of standard MD. As examples, we performed two simulations of a million-atom system to model the diffusion and clustering of Pt adatoms on a large patch of the Pt(100) surface for 80 μs and 160 µs.