Enhanced sampling algorithms

Alchemical Enhanced Sampling with Optimized Phase Space Overlap

An alchemical enhanced sampling (ACES) method has recently been introduced to facilitate importance sampling in free energy simulations. The method achieves enhanced sampling from Hamiltonian replica exchange within a dual topology framework while utilizing new smoothstep softcore potentials. A common sampling problem encountered in lead optimization is the functionalization of aromatic rings that exhibit distinct conformational preferences when interacting with the protein. It is difficult to converge the distribution of ring conformations due to the long time scale of ring flipping events; however, the ACES method addresses this issue by modeling the syn and anti ring conformations within a dual topology. ACES thereby samples the conformer distributions by alchemically tunneling between states, as opposed to traversing a physical pathway with a high rotational barrier. We demonstrate the use of ACES to overcome conformational sampling issues involving ring flipping in ML300-derived noncovalent inhibitors of SARS-CoV-2 Main Protease (Mpro). The demonstrations explore how the use of replica exchange and the choice of softcore selection affects the convergence of the ring conformation distributions. Furthermore, we examine how the accuracy of the calculated free energies is affected by the degree of phase space overlap (PSO) between adjacent states (i.e., between neighboring λ-windows) and the Hamiltonian replica exchange acceptance ratios. Both of these factors are sensitive to the spacing between the intermediate states. We introduce a new method for choosing a schedule of λ values. The method analyzes short "burn-in" simulations to construct a 2D map of the nonlocal PSO. The schedule is obtained by optimizing an alchemical pathway on the 2D map that equalizes the PSO between the λ intervals. The optimized phase space overlap λ-spacing method (Opt-PSO) leads to more numerous end-to-end single passes and round trips due to the correlation between PSO and Hamiltonian replica exchange acceptance ratios. The improved exchange statistics enhance the efficiency of ACES method. The method has been implemented into the FE-ToolKit software package, which is freely available.

Molecular Dynamics and Other HPC Simulations for Drug Discovery

High performance computing (HPC) is taking an increasingly important place in drug discovery. It makes possible the simulation of complex biochemical systems with high precision in a short time, thanks to the use of sophisticated algorithms. It promotes the advancement of knowledge in fields that are inaccessible or difficult to access through experimentation and it contributes to accelerating the discovery of drugs for unmet medical needs while reducing costs. Herein, we report how computational performance has evolved over the past years, and then we detail three domains where HPC is essential. Molecular dynamics (MD) is commonly used to explore the flexibility of proteins, thus generating a better understanding of different possible approaches to modulate their activity. Modeling and simulation of biopolymer complexes enables the study of protein-protein interactions (PPI) in healthy and disease states, thus helping the identification of targets of pharmacological interest. Virtual screening (VS) also benefits from HPC to predict in a short time, among millions or billions of virtual chemical compounds, the best potential ligands that will be tested in relevant assays to start a rational drug design process.

When MELD Meets GaMD: Accelerating Biomolecular Landscape Exploration

We introduce Gaussian accelerated MELD (GaMELD) as a new method for exploring the energy landscape of biomolecules. GaMELD combines the strengths of Gaussian accelerated molecular dynamics (GaMD) and modeling employing limited data (MELD) to navigate complex energy landscapes. MELD uses replica-exchange molecular simulations to integrate limited and uncertain data into simulations via Bayesian inference. MELD has been successfully applied to problems of structure prediction like protein folding and complex structure prediction. However, the computational cost for MELD simulations has limited its broader applicability. The synergy of GaMD and MELD surmounts this limitation efficiently sampling the energy landscape at a lower computational cost (reducing the computational cost by a factor of 2 to six). Effectively, GaMD is used to shift energy distributions along replicas to increase the overlap in energy distributions across replicas, facilitating a random walk in replica space. We tested GaMELD on a benchmark set of 12 small proteins that have been previously studied through MELD and conventional MD. GaMELD consistently achieves accurate predictions with fewer replicas. By increasing the efficacy of replica exchange, GaMELD effectively accelerates convergence in the conformational space, enabling improved sampling across a diverse set of systems.

Modern Alchemical Free Energy Methods for Drug Discovery Explained

This Perspective provides a contextual explanation of the current state-of-the-art alchemical free energy methods and their role in drug discovery as well as highlights select emerging technologies. The narrative attempts to answer basic questions about what goes on "under the hood" in free energy simulations and provide general guidelines for how to run simulations and analyze the results. It is the hope that this work will provide a valuable introduction to students and scientists in the field.

A New Approach for Estimating the Free Energy Differences among Multiple Thermodynamic States in Statistical Simulations

In this letter, a new approach to compute free energy differences (FEDs) between multiple thermodynamics states is introduced. The method directly uses energy probability densities, which can be extracted with high accuracy from equilibrium simulations to obtain FEDs. Methods in current use, such as Bennett acceptance ratio (BAR), its multistate generalization (MBAR), or the weighted histogram analysis method (WHAM), require iterative solution of nonlinear equations which are known to be slowly convergent. The equations providing MBAR FEDs are identical to those derived earlier by Souaille and Roux in a method that has become known informally as "binless WHAM". In contrast, we obtain FEDs by solution of linear equations. For the classic two-state problem, the statistical error of our method, solving linear equations, is shown analytically to match that of BAR under common conditions.

ACES: Optimized Alchemically Enhanced Sampling

We present an alchemical enhanced sampling (ACES) method implemented in the GPU-accelerated AMBER free energy MD engine. The methods hinges on the creation of an "enhanced sampling state" by reducing or eliminating selected potential energy terms and interactions that lead to kinetic traps and conformational barriers while maintaining those terms that curtail the need to otherwise sample large volumes of phase space. For example, the enhanced sampling state might involve transforming regions of a ligand and/or protein side chain into a noninteracting "dummy state" with internal electrostatics and torsion angle terms turned off. The enhanced sampling state is connected to a real-state end point through a Hamiltonian replica exchange (HREMD) framework that is facilitated by newly developed alchemical transformation pathways and smoothstep softcore potentials. This creates a counterdiffusion of real and enhanced-sampling states along the HREMD network. The effect of a differential response of the environment to the real and enhanced-sampling states is minimized by leveraging the dual topology framework in AMBER to construct a counterbalancing HREMD network in the opposite alchemical direction with the same (or similar) real and enhanced sampling states at inverted end points. The method has been demonstrated in a series of test cases of increasing complexity where traditional MD, and in several cases alternative REST2-like enhanced sampling methods, are shown to fail. The hydration free energy for acetic acid was shown to be independent of the starting conformation, and the values for four additional edge case molecules from the FreeSolv database were shown to have a significantly closer agreement with experiment using ACES. The method was further able to handle different rotamer states in a Cdk2 ligand identified as fractionally occupied in crystal structures. Finally, ACES was applied to T4-lysozyme and demonstrated that the side chain distribution of V111χ1 could be reliably reproduced for the apo state, bound to p-xylene, and in p-xylene→ benzene transformations. In these cases, the ACES method is shown to be highly robust and superior to a REST2-like enhanced sampling implementation alone.

Temporally Coherent Backmapping of Molecular Trajectories From Coarse-Grained to Atomistic Resolution

Coarse-graining offers a means to extend the achievable time and length scales of molecular dynamics simulations beyond what is practically possible in the atomistic regime. Sampling molecular configurations of interest can be done efficiently using coarse-grained simulations, from which meaningful physicochemical information can be inferred if the corresponding all-atom configurations are reconstructed. However, this procedure of backmapping to reintroduce the lost atomistic detail into coarse-grain structures has proven a challenging task due to the many feasible atomistic configurations that can be associated with one coarse-grain structure. Existing backmapping methods are strictly frame-based, relying on either heuristics to replace coarse-grain particles with atomic fragments and subsequent relaxation or parametrized models to propose atomic coordinates separately and independently for each coarse-grain structure. These approaches neglect information from previous trajectory frames that is critical to ensuring temporal coherence of the backmapped trajectory, while also offering information potentially helpful to producing higher-fidelity atomic reconstructions. In this work, we present a deep learning-enabled data-driven approach for temporally coherent backmapping that explicitly incorporates information from preceding trajectory structures. Our method trains a conditional variational autoencoder to nondeterministically reconstruct atomistic detail conditioned on both the target coarse-grain configuration and the previously reconstructed atomistic configuration. We demonstrate our backmapping approach on two exemplar biomolecular systems: alanine dipeptide and the miniprotein chignolin. We show that our backmapped trajectories accurately recover the structural, thermodynamic, and kinetic properties of the atomistic trajectory data.

Unraveling the Role of Charge Patterning in the Micellar Structure of Sequence-Defined Amphiphilic Peptoid Oligomers by Molecular Dynamics Simulations

Electrostatic interactions play a significant role in regulating biological systems and have received increasing attention due to their usefulness in designing advanced stimulus-responsive materials. Polypeptoids are highly tunable N-substituted peptidomimetic polymers that lack backbone hydrogen bonding and chirality. Therefore, polypeptoids are suitable systems to study the effect of noncovalent interactions of substituents without complications of backbone intramolecular and intermolecular hydrogen bonding. In this study, all-atom molecular dynamics (MD) simulations were performed on micelles formed by a series of sequence-defined ionic polypeptoid block copolymers consisting of a hydrophobic segment and a hydrophilic segment in an aqueous solution. By combining the results from MD simulations and experimental small-angle neutron scattering data, further insights were gained into the internal structure of the formed polypeptoid micelles, which is not always directly accessible from experiments. In addition, information was gained into the physics of the noncovalent interactions responsible for the self-assembly of weakly charged polypeptoids in an aqueous solution. While the aggregation number is governed by electrostatic repulsion of the negatively charged carboxylate (COO^–) substituents on the polypeptoid chain within the micelle, MD simulations indicate that the position of the charge on singly charged chains mediates the shape of the micelle through the charge–dipole interactions between the COO^– substituent and the surrounding water. Therefore, the polypeptoid micelles formed from the single-charged series offer the possibility for tailorable micelle shapes. In contrast, the polypeptoid micelles formed from the triple-charged series are characterized by more pronounced electrostatic repulsion that competes with more significant charge–sodium interactions, making it difficult to predict the shape of the micelles. This work has helped further develop design principles for the shape and structure of self-assembled micelles by controlling the position of charged moieties on the backbone of polypeptoid block copolymers.

Effect of Lithium Drug on Binding Affinities of Glycogen Synthase Kinase-3 β to Its Network Partners: A New Computational Approach

Finding new methods to study the effect of small molecules on protein interaction networks provides us with invaluable tools in the fields of pharmacodynamics and drug design. Lithium is an antimanic drug that has been used for the treatment of bipolar disorder for more than 60 years. Here, we utilized a new approach to study the effect of lithium as a drug on the protein interaction network of GSK-3β as a hub protein and computed the affinities of GSK-3β to its partners in the presence of lithium or sodium ions. For this purpose, ensembles of GSK-3β protein structures were created in the presence of either lithium or sodium ions using adaptive tempering molecular dynamics simulations. The protein binding patches of GSK-3β for its partners were determined, and finally, the affinity of each binding patch to the related partner was computed for structures of ensembles using a monomer-based approach. Besides, by comparing structural dynamics of GSK-3β during MD simulations in the presence of LiCl and NaCl, we suggested a new mechanism for the inhibitory effect of lithium on GSK-3β.

To the fast calculation of the solvation free energy. Combining expanded ensembles with L2MC

A more efficient version of the Expanded Ensembles method for calculation of free energy in molecular-mechanical simulations is proposed. The method is based on the Horowitz L2MC approach to accelerate movement along the "alchemical" coordinate. It is possible to achieve the same efficiency of the algorithm both with the optimal number of "windows" and with a larger number of them compared to the original algorithm. Since the optimal number of windows is unknown a priory, the proposed algorithm is more robust than the traditional one. We can choose the number of windows in excess and do not worry about the loss of efficiency. We illustrate the method's efficiency with the computation of the hydration free energies of pyridine and water.

Enhanced solvation force extrapolation for speeding up molecular dynamics simulations of complex biochemical liquids

We propose an enhanced approach to the extrapolation of mean potential forces acting on atoms of solute macromolecules due to their interactions with solvent atoms in complex biochemical liquids. It improves and extends our previous extrapolation schemes by additionally including new techniques such as an exponential scaling transformation of coordinate space with weights complemented by an automatically adjusted balancing between the least square minimization of force deviations and the norm of expansion coefficients in the approximation. The expensive mean potential forces are treated in terms of the three-dimensional reference interaction site model with Kovalenko-Hirata closure molecular theory of solvation. During the dynamics, they are calculated only after every long (outer) time interval, i.e., quite rarely to reduce the computational costs. At much shorter (inner) time steps, these forces are extrapolated on the basis of their outer values. The equations of motion are then solved using a multiple time step integration within an optimized isokinetic Nosé-Hoover chain thermostat. The new approach is applied to molecular dynamics simulations of various systems consisting of solvated organic and biomolecules of different complexity. For example, we consider hydrated alanine dipeptide, asphaltene in toluene solvent, miniprotein 1L2Y, and protein G in aqueous solution. It is shown that in all these cases, the enhanced extrapolation provides much better accuracy of the solvation force approximation than the existing approaches. As a result, it can be used with much larger outer time steps, leading to a significant speedup of the simulations.

Determination of the structural ensemble of the molten globule state of a protein by computer simulations

We have used computer simulations to investigate the structural nature of the molten globule (MG) state of canine milk lysozyme. To sample the conformational space efficiently, we performed replica-exchange umbrella sampling simulations with the radius of gyration as a reaction coordinate. We applied the Weighted Histogram Analysis Method to the trajectory of the simulations to obtain the potential of mean force, from which we identified representative structures corresponding to local minima in the free energy surface. The representative structures obtained in this way are in accord with the characteristics of the MG state reported previously by experimental studies. We conjecture that the MG state comprises a series of partially structured states undergoing relatively fast conformational interchange.

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.

Quantum chemical replica-exchange umbrella sampling molecular dynamics simulations reveal the formation mechanism of iron phthalocyanine from iron and phthalonitrile

Phthalocyanine (Pc) and its metal complexes (MPcs) have been used industrially since their discovery in the early 20th century. The phthalonitrile (PN) method is a well-known synthesis method in which Pc or MPc can be afforded by heating a mixture of PN and metal powders over 280 °C with only moderate yield. However, the formation mechanism of the phthalocyanines and the intermediate stages of this seemingly simple reaction have yet to be fully understood. To study this mechanism computationally, we carried out quantum chemical molecular dynamics (MD) simulations based on the density-functional tight-binding (DFTB) method, applying the replica-exchange umbrella sampling (REUS) method, starting from four PN molecules and one iron atom. The DFTB-REUS-MD simulations successfully yielded FePc, and a metastable structure very similar to FePc but with a reactive nitrene unit was also identified that might explain the incomplete conversion of the reactants into FePc. Analysis of the MD trajectories reveals a three-step FePc formation mechanism for the PN method.

Towards a coarse-grained model of the peptoid backbone: the case of N,N-dimethylacetamide

Coarse-grained model of DMA, containing the basic motif of the peptoid backbone, based on short ranged many-body ranged interactions.

Replica Exchange Molecular Dynamics: A Practical Application Protocol with Solutions to Common Problems and a Peptide Aggregation and Self-Assembly Example

Protein aggregation is associated with many human diseases such as Alzheimer's disease (AD), Parkinson's disease (PD), and type II diabetes (T2D). Understanding the molecular mechanism of protein aggregation is essential for therapy development. Molecular dynamics (MD) simulations have been shown as powerful tools to study protein aggregation. However, conventional MD simulations can hardly sample the whole conformational space of complex protein systems within acceptable simulation time as it can be easily trapped in local minimum-energy states. Many enhanced sampling methods have been developed. Among these, the replica exchange molecular dynamics (REMD) method has gained great popularity. By combining MD simulation with the Monte Carlo algorithm, the REMD method is capable of overcoming high energy-barriers easily and of sampling sufficiently the conformational space of proteins. In this chapter, we present a brief introduction to REMD method and a practical application protocol with a case study of the dimerization of the 11-25 fragment of human islet amyloid polypeptide (hIAPP(11-25)), using the GROMACS software. We also provide solutions to problems that are often encountered in practical use, and provide some useful scripts/commands from our research that can be easily adapted to other systems.

Approaching Pharmacological Space: Events and Components

With a view to introducing the concept of pharmacological space and its potential applications in investigating and predicting the toxic mechanisms of xenobiotics, this opening chapter describes the logical relations between conformational behavior, physicochemical properties and binding spaces, which are seen as the three key elements composing the pharmacological space. While the concept of conformational space is routinely used to encode molecular flexibility, the concepts of property spaces and, particularly, of binding spaces are more innovative. Indeed, their descriptors can find fruitful applications (a) in describing the dynamic adaptability a given ligand experiences when inserted into a specific environment, and (b) in parameterizing the flexibility a ligand retains when bound to a biological target. Overall, these descriptors can conveniently account for the often disregarded entropic factors and as such they prove successful when inserted in ligand- or structure-based predictive models. Notably, and although binding space parameters can clearly be derived from MD simulations, the chapter will illustrate how docking calculations, despite their static nature, are able to evaluate ligand's flexibility by analyzing several poses for each ligand. Such an approach, which represents the founding core of the binding space concept, can find various applications in which the related descriptors show an impressive enhancing effect on the statistical performances of the resulting predictive models.

SAXS-Oriented Ensemble Refinement of Flexible Biomolecules

The conformational flexibility of a biomolecule may play a crucial role in its biological function. Small-angle x-ray scattering (SAXS) is a very popular technique for characterizing biomolecule flexibility. It can be used to determine a possible structural ensemble of the biomolecule in solution with the aid of a computer simulation. In this article, we present a tool written in Python, which iteratively runs multiple independent enhanced sampling simulations such as amplified collective motions and accelerated molecular dynamics, and an ensemble optimization method to drive the biomolecule toward an ensemble that fits the SAXS data well. The tool has been validated with a protein and an RNA system, i.e., the tandem WW domains of formin-binding protein 21 and the aptamer domain of SAM-1 riboswitch, respectively. These Python scripts are user-friendly and can be easily modified if a different simulation engine is preferred.

Biasing Simulations of DNA Base Pair Parameters with Application to Propellor Twisting in AT/AT, AA/TT, and AC/GT Steps and Their Uracil Analogs

An accurate and efficient implementation of the six DNA base pair parameters as order parameters for enhanced sampling simulations is presented. The parameter definitions are defined by vector algebra operations on a reduced atomic set of the base pair, and correlate very well with standard definitions. Application of the model is illustrated by umbrella sampling simulations of propeller twisting within AT/AT, AA/TT, and AC/GT steps and their uracil analogs. Strong correlations are found between propeller twisting and a number of conformational parameters, including buckle, opening, BI/BII backbone configuration, and sugar puckering. The thymine methyl group is observed to notably alter the local conformational free energy landscape, with effects within and directly upstream of the thymine containing base pair.

Enhanced Sampling of an Atomic Model with Hybrid Nonequilibrium Molecular Dynamics-Monte Carlo Simulations Guided by a Coarse-Grained Model

Molecular dynamics (MD) trajectories based on a classical equation of motion provide a straightforward, albeit somewhat inefficient approach, to explore and sample the configurational space of a complex molecular system. While a broad range of techniques can be used to accelerate and enhance the sampling efficiency of classical simulations, only algorithms that are consistent with the Boltzmann equilibrium distribution yield a proper statistical mechanical computational framework. Here, a multiscale hybrid algorithm relying simultaneously on all-atom fine-grained (FG) and coarse-grained (CG) representations of a system is designed to improve sampling efficiency by combining the strength of nonequilibrium molecular dynamics (neMD) and Metropolis Monte Carlo (MC). This CG-guided hybrid neMD-MC algorithm comprises six steps: (1) a FG configuration of an atomic system is dynamically propagated for some period of time using equilibrium MD; (2) the resulting FG configuration is mapped onto a simplified CG model; (3) the CG model is propagated for a brief time interval to yield a new CG configuration; (4) the resulting CG configuration is used as a target to guide the evolution of the FG system; (5) the FG configuration (from step 1) is driven via a nonequilibrium MD (neMD) simulation toward the CG target; (6) the resulting FG configuration at the end of the neMD trajectory is then accepted or rejected according to a Metropolis criterion before returning to step 1. A symmetric two-ends momentum reversal prescription is used for the neMD trajectories of the FG system to guarantee that the CG-guided hybrid neMD-MC algorithm obeys microscopic detailed balance and rigorously yields the equilibrium Boltzmann distribution. The enhanced sampling achieved with the method is illustrated with a model system with hindered diffusion and explicit-solvent peptide simulations. Illustrative tests indicate that the method can yield a speedup of about 80 times for the model system and up to 21 times for polyalanine and (AAQAA)₃ in water.

The mechanism of water/ion exchange at a protein surface: a weakly bound chloride in Helicobacter pylori apoflavodoxin

Binding/unbinding of small ligands, such as ions, to/from proteins influences biochemical processes such as protein folding, enzyme catalysis or protein/ligand recognition. We have investigated the mechanism of chloride/water exchange at a protein surface (that of the apoflavodoxin from Helicobacter pylori) using classical all-atom molecular dynamics simulations. They reveal a variety of chloride exit routes and residence times; the latter is related to specific coordination modes of the anion. The role of solvent molecules in the mechanism of chloride unbinding has been studied in detail. We see no temporary increase in chloride coordination along the release process. Instead, the coordination of new water molecules takes place in most cases after the chloride/protein atom release event has begun. Moreover, the distribution function of water entrance events into the first chloride solvation shell peaks after chloride protein atom dissociation events. All these observations together seem to indicate that water molecules simply fill the vacancies left by the previously coordinating protein residues. We thus propose a step-by-step dissociation pathway in which protein/chloride interactions gradually break down before new water molecules progressively fill the vacant positions left by protein atoms. As observed for other systems, water molecules associated with bound chloride or with protein atoms have longer residence times than those bound to the free anion. The implications of the exchange mechanism proposed for the binding of the FMN (Flavin Mononucleotide) protein cofactor are discussed.

Anisotropy of B-DNA Groove Bending

DNA bending is critical for DNA packaging, recognition, and repair, and occurs toward either the major or the minor groove. The anisotropy of B-DNA groove bending was quantified for eight DNA sequences by free energy simulations employing a novel reaction coordinate. The simulations show that bending toward the major groove is preferred for non-A-tracts while the A-tract has a high tendency of bending toward the minor groove. Persistence lengths were generally larger for bending toward the minor groove, which is thought to originate from differences in groove hydration. While this difference in stiffness is one of the factors determining the overall preference of bending direction, the dominant contribution is shown to be a free energy offset between major and minor groove bending. The data suggests that, for the A-tract, this offset is largely determined by inherent structural properties, while differences in groove hydration play a large role for non-A-tracts. By quantifying the energetics of DNA groove bending and rationalizing the origins of the anisotropy, the calculations provide important new insights into a key biological process.

More than a look into a crystal ball: protein structure elucidation guided by molecular dynamics simulations

The 'form follows function' principle implies that a structural determination of protein structures is indispensable to understand proteins in their biological roles. However, experimental methods still show shortcomings in the description of the dynamic properties of proteins. Therefore, molecular dynamics (MD) simulations represent an essential tool for structural biology to investigate proteins as flexible and dynamic entities. Here, we will give an overview on the impact of MD simulations on structural investigations, including studies that aim at a prediction of protein-folding pathways, protein-assembly processes and the sampling of conformational space by computational means.

A critical assessment of hidden markov model sub-optimal sampling strategies applied to the generation of peptide 3D models

Hidden Markov Model derived structural alphabets are a probabilistic framework in which the complete conformational space of a peptidic chain is described in terms of probability distributions that can be sampled to identify conformations of largest probabilities. Here, we assess how three strategies to sample sub-optimal conformations-Viterbi k-best, forward backtrack and a taboo sampling approach-can lead to the efficient generation of peptide conformations. We show that the diversity of sampling is essential to compensate biases introduced in the estimates of the probabilities, and we find that only the forward backtrack and a taboo sampling strategies can efficiently generate native or near-native models. Finally, we also find such approaches are as efficient as former protocols, while being one order of magnitude faster, opening the door to the large scale de novo modeling of peptides and mini-proteins. © 2016 Wiley Periodicals, Inc.

Characterization of the three-dimensional free energy manifold for the uracil ribonucleoside from asynchronous replica exchange simulations

Replica exchange molecular dynamics has emerged as a powerful tool for efficiently sampling free energy landscapes for conformational and chemical transitions. However, daunting challenges remain in efficiently getting such simulations to scale to the very large number of replicas required to address problems in state spaces beyond two dimensions. The development of enabling technology to carry out such simulations is in its infancy, and thus it remains an open question as to which applications demand extension into higher dimensions. In the present work, we explore this problem space by applying asynchronous Hamiltonian replica exchange molecular dynamics with a combined quantum mechanical/molecular mechanical potential to explore the conformational space for a simple ribonucleoside. This is done using a newly developed software framework capable of executing >3,000 replicas with only enough resources to run 2,000 simultaneously. This may not be possible with traditional synchronous replica exchange approaches. Our results demonstrate 1.) the necessity of high dimensional sampling simulations for biological systems, even as simple as a single ribonucleoside, and 2.) the utility of asynchronous exchange protocols in managing simultaneous resource requirements expected in high dimensional sampling simulations. It is expected that more complicated systems will only increase in computational demand and complexity, and thus the reported asynchronous approach may be increasingly beneficial in order to make such applications available to a broad range of computational scientists.

Dynamical traps in Wang-Landau sampling of continuous systems: Mechanism and solution

We study the mechanism behind dynamical trappings experienced during Wang-Landau sampling of continuous systems reported by several authors. Trapping is caused by the random walker coming close to a local energy extremum, although the mechanism is different from that of the critical slowing-down encountered in conventional molecular dynamics or Monte Carlo simulations. When trapped, the random walker misses the entire or even several stages of Wang-Landau modification factor reduction, leading to inadequate sampling of the configuration space and a rough density of states, even though the modification factor has been reduced to very small values. Trapping is dependent on specific systems, the choice of energy bins, and the Monte Carlo step size, making it highly unpredictable. A general, simple, and effective solution is proposed where the configurations of multiple parallel Wang-Landau trajectories are interswapped to prevent trapping. We also explain why swapping frees the random walker from such traps. The efficacy of the proposed algorithm is demonstrated.

Liquid-solid and solid-solid phase transition of monolayer water: high-density rhombic monolayer ice

Liquid-solid and solid-solid phase transitions of a monolayer water confined between two parallel hydrophobic surfaces are studied by molecular dynamics simulations. The solid phase considered is the high-density rhombic monolayer ice. Based on the computed free energy surface, it is found that at a certain width of the slit nanopore, the monolayer water exhibits not only a high freezing point but also a low energy barrier to crystallization. Moreover, through analyzing the oxygen-hydrogen-oxygen angle distribution and oxygen-hydrogen radial distribution, the high-density monolayer ice is classified as either a flat ice or a puckered ice. The transition between a flat ice and a puckered ice reflects a trade-off between the water-wall interactions and the electrostatic interactions among water molecules.

Feedback from network states generates variability in a probabilistic olfactory circuit

Variability is a prominent feature of behavior and is an active element of certain behavioral strategies. To understand how neuronal circuits control variability, we examined the propagation of sensory information in a chemotaxis circuit of C. elegans where discrete sensory inputs can drive a probabilistic behavioral response. Olfactory neurons respond to odor stimuli with rapid and reliable changes in activity, but downstream AIB interneurons respond with a probabilistic delay. The interneuron response to odor depends on the collective activity of multiple neurons-AIB, RIM, and AVA-when the odor stimulus arrives. Certain activity states of the network correlate with reliable responses to odor stimuli. Artificially generating these activity states by modifying neuronal activity increases the reliability of odor responses in interneurons and the reliability of the behavioral response to odor. The integration of sensory information with network states may represent a general mechanism for generating variability in behavior.

Investigating drug-target association and dissociation mechanisms using metadynamics-based algorithms

CONSPECTUS: This Account highlights recent advances and discusses major challenges in the field of drug-target recognition, binding, and unbinding studied using metadynamics-based approaches, with particular emphasis on their role in structure-based design. Computational chemistry has significantly contributed to drug design and optimization in an extremely broad range of areas, including prediction of target druggability and drug likeness, de novo design, fragment screening, ligand docking, estimation of binding affinity, and modulation of ADMET (absorption, distribution, metabolism, excretion, toxicity) properties. Computationally driven drug discovery must continuously adapt to keep pace with the evolving knowledge of the factors that modulate the pharmacological action of drugs. There is thus an urgent need for novel computational approaches that integrate the vast amount of complex information currently available for small (bio)organic compounds, biologically relevant targets and their complexes, while also accounting accurately for the thermodynamics and kinetics of drug-target association, the intrinsic dynamical behavior of biomolecular systems, and the complexity of protein-protein networks. Understanding the mechanism of drug binding to and unbinding from biological targets is fundamental for optimizing lead compounds and designing novel biologically active ones. One major challenge is the accurate description of the conformational complexity prior to and upon formation of drug-target complexes. Recently, enhanced sampling methods, including metadynamics and related approaches, have been successfully applied to investigate complex mechanisms of drugs binding to flexible targets. Metadynamics is a family of enhanced sampling techniques aimed at enhancing the rare events and reconstructing the underlying free energy landscape as a function of a set of order parameters, usually referred to as collective variables. Studies of drug binding mechanisms have predicted the most probable association and dissociation pathways and the related binding free energy profile. In addition, the availability of an efficient open-source implementation, running on cost-effective GPU (i.e., graphical processor unit) architectures, has considerably decreased the learning curve and the computational costs of the methods, and increased their adoption by the community. Here, we review the recent contributions of metadynamics and other enhanced sampling methods to the field of drug-target recognition and binding. We discuss how metadynamics has been used to search for transition states, to predict binding and unbinding paths, to treat conformational flexibility, and to compute free energy profiles. We highlight the importance of such predictions in drug discovery. Major challenges in the field and possible solutions will finally be discussed.

Computational Recipe for Efficient Description of Large-Scale Conformational Changes in Biomolecular Systems

Characterizing large-scale structural transitions in biomolecular systems poses major technical challenges to both experimental and computational approaches. On the computational side, efficient sampling of the configuration space along the transition pathway remains the most daunting challenge. Recognizing this issue, we introduce a knowledge-based computational approach toward describing large-scale conformational transitions using (i) nonequilibrium, driven simulations combined with work measurements and (ii) free energy calculations using empirically optimized biasing protocols. The first part is based on designing mechanistically relevant, system-specific reaction coordinates whose usefulness and applicability in inducing the transition of interest are examined using knowledge-based, qualitative assessments along with nonequilirbrium work measurements which provide an empirical framework for optimizing the biasing protocol. The second part employs the optimized biasing protocol resulting from the first part to initiate free energy calculations and characterize the transition quantitatively. Using a biasing protocol fine-tuned to a particular transition not only improves the accuracy of the resulting free energies but also speeds up the convergence. The efficiency of the sampling will be assessed by employing dimensionality reduction techniques to help detect possible flaws and provide potential improvements in the design of the biasing protocol. Structural transition of a membrane transporter will be used as an example to illustrate the workings of the proposed approach.

Statistical mechanics of the denatured state of a protein using replica-averaged metadynamics

The characterization of denatured states of proteins is challenging because the lack of permanent structure in these states makes it difficult to apply to them standard methods of structural biology. In this work we use all-atom replica-averaged metadynamics (RAM) simulations with NMR chemical shift restraints to determine an ensemble of structures representing an acid-denatured state of the 86-residue protein ACBP. This approach has enabled us to reach convergence in the free energy landscape calculations, obtaining an ensemble of structures in relatively accurate agreement with independent experimental data used for validation. By observing at atomistic resolution the transient formation of native and non-native structures in this acid-denatured state of ACBP, we rationalize the effects of single-point mutations on the folding rate, stability, and transition-state structures of this protein, thus characterizing the role of the unfolded state in determining the folding process.

Recognition of the HIV capsid by the TRIM5α restriction factor is mediated by a subset of pre-existing conformations of the TRIM5α SPRY domain

The binding of the TRIM5α restriction factor to the HIV capsid is mediated by the C-terminal SPRY domain of TRIM5α. Atomic-level details of this host-pathogen interaction, which involves mobile variable loops of the SPRY domain, remain unclear. Some of the key determinants of restriction are encompassed by the long and disordered v1 loop of the SPRY domain. We applied molecular modeling to elucidate the conformational repertoire of the v1 loop and its role in the interaction with the capsid. All-atom replica exchange molecular dynamics revealed multiple transient, interconverting states of the v1 loop consistent with the intrinsic disorder observed experimentally. The docking of the SPRY conformations representing 10 most populated states onto the high-resolution model of the assembled HIV-1 capsid revealed that a subset of v1 conformations produced plausible binding poses, in which the SPRY domain binds close to the pseudo-2-fold symmetry axis and the v1 loop spans the interhexamer gap. Such binding mode is well supported by the NMR binding data and known escape mutants. We speculate that the binding mode that involves interaction of the capsid with a subset of preexisting SPRY conformations arising from the intrinsic disorder of the v1 loop may explain the remarkable ability of TRIM5α to resist viral evasion by mutagenesis and to restrict divergent retroviruses.

Accelerated Molecular Dynamics Simulations with the AMOEBA Polarizable Force Field on Graphics Processing Units

The accelerated molecular dynamics (aMD) method has recently been shown to enhance the sampling of biomolecules in molecular dynamics (MD) simulations, often by several orders of magnitude. Here, we describe an implementation of the aMD method for the OpenMM application layer that takes full advantage of graphics processing units (GPUs) computing. The aMD method is shown to work in combination with the AMOEBA polarizable force field (AMOEBA-aMD), allowing the simulation of long time-scale events with a polarizable force field. Benchmarks are provided to show that the AMOEBA-aMD method is efficiently implemented and produces accurate results in its standard parametrization. For the BPTI protein, we demonstrate that the protein structure described with AMOEBA remains stable even on the extended time scales accessed at high levels of accelerations. For the DNA repair metalloenzyme endonuclease IV, we show that the use of the AMOEBA force field is a significant improvement over fixed charged models for describing the enzyme active-site. The new AMOEBA-aMD method is publicly available (http://wiki.simtk.org/openmm/VirtualRepository) and promises to be interesting for studying complex systems that can benefit from both the use of a polarizable force field and enhanced sampling.

Pressure-Induced Helical Structure of a Peptide Studied by Simulated Tempering Molecular Dynamics Simulations

It is known experimentally that an AK16 peptide forms more α-helix structures with increasing pressure while proteins unfold in general. In order to understand this abnormality, molecular dynamics (MD) simulations with the simulated tempering method for the isobaric-isothermal ensemble were performed in a wide pressure range from 1.0 × 10(-4) GPa to 1.4 GPa. From the results of the simulations, it is found that the fraction of the folded state decreases once and increases after that with increasing pressure. The partial molar volume change from the folded state to unfolded state increases monotonically from a negative value to a positive value with pressure. The behavior under high pressure conditions is consistent with the experimental results. The radius of gyration of highly helical structures decreases with increasing pressure, which indicates that the helix structure shrinks with pressure. This is the reason why the fraction of the folded state increases as pressure increases.