Achieving ergodic sampling using replica-exchange free-energy calculations

Predictions from First-Principles of Membrane Permeability to Small Molecules: How Useful Are They in Practice?

Predicting from first-principles the rate of passive permeation of small molecules across the biological membrane represents a promising strategy for screening lead compounds upstream in the drug-discovery and development pipeline. One popular avenue for the estimation of permeation rates rests on computer simulations in conjunction with the inhomogeneous solubility-diffusion model, which requires the determination of the free-energy change and position-dependent diffusivity of the substrate along the translocation pathway through the lipid bilayer. In this Perspective, we will clarify the physical meaning of the membrane permeability inferred from such computer simulations, and how theoretical predictions actually relate to what is commonly measured experimentally. We will also examine why these calculations remain both technically challenging and overly computationally expensive, which has hitherto precluded their routine use in nonacademic settings. We finally synopsize possible research directions to meet these challenges, increase the predictive power of physics-based rates of passive permeation, and, by ricochet, improve their practical usefulness.

Enhancing sampling with free-energy calculations

In recent years, considerable progress has been made to enhance sampling and help address biological questions, including, but not limited to conformational transitions in biomolecules and protein-ligand reversible binding, hitherto intractable by brute-force computer simulations. Many of these advances result from the development of a palette of methods aimed at exploring rare events through reliable free-energy calculations. The advent of new, often conceptually related methods has also rendered difficult the choice of the best suited option for a given problem. Here, we focus on geometrical transformations and algorithms designed to enhance sampling along adequately chosen progress variables, tracing their theoretical foundations, and showing how they are connected and can be blended together for improved performance.

Overcoming Free-Energy Barriers with a Seamless Combination of a Biasing Force and a Collective Variable-Independent Boost Potential

Amid collective-variable (CV)-based importance-sampling algorithms, a hybrid of the extended adaptive biasing force and the well-tempered metadynamics algorithms (WTM-eABF) has proven particularly cost-effective for exploring the rugged free-energy landscapes that underlie biological processes. However, as an inherently CV-based algorithm, this hybrid scheme does not explicitly accelerate sampling in the space orthogonal to the chosen CVs, thereby limiting its efficiency and accuracy, most notably in those cases where the slow degrees of freedom of the process at hand are not accounted for in the model transition coordinate. Here, inspired by Gaussian-accelerated molecular dynamics (GaMD), we introduce the same CV-independent harmonic boost potential into WTM-eABF, yielding a hybrid algorithm coined GaWTM-eABF. This algorithm leans on WTM-eABF to explore the transition coordinate with a GaMD-mollified potential and recovers the unbiased free-energy landscape through thermodynamic integration followed by proper reweighting. As illustrated in our numerical tests, GaWTM-eABF effectively overcomes the free-energy barriers in orthogonal space and correctly recovers the unbiased potential of mean force (PMF). Furthermore, applying both GaWTM-eABF and WTM-eABF to two biologically relevant processes, namely, the reversible folding of (i) deca-alanine and (ii) chignolin, our results indicate that GaWTM-eABF reduces the uncertainty in the PMF calculation and converges appreciably faster than WTM-eABF. Obviating the need of multiple-copy strategies, GaWTM-eABF is a robust, computationally efficient algorithm to surmount the free-energy barriers in orthogonal space and maps with utmost fidelity the free-energy landscape along selections of CVs. Moreover, our strategy that combines WTM-eABF with GaMD can be easily extended to other biasing-force algorithms.

Atomistic Basis of Microtubule Dynamic Instability Assessed Via Multiscale Modeling

Microtubule “dynamic instability,” the abrupt switching from assembly to disassembly caused by the hydrolysis of GTP to GDP within the β subunit of the αβ-tubulin heterodimer, is necessary for vital cellular processes such as mitosis and migration. Despite existing high-resolution structural data, the key mechanochemical differences between the GTP and GDP states that mediate dynamic instability behavior remain unclear. Starting with a published atomic-level structure as an input, we used multiscale modeling to find that GTP hydrolysis results in both longitudinal bond weakening (~ 4 k_BT) and an outward bending preference (~ 1.5 k_BT) to both drive dynamic instability and give rise to the microtubule tip structures previously observed by light and electron microscopy. More generally, our study provides an example where atomic level structural information is used as the sole input to predict cellular level dynamics without parameter adjustment.

Potential of Matrix Metalloproteinase Inhibitors for the Treatment of Local Tissue Damage Induced by a Type P-I Snake Venom Metalloproteinase

Snake bite envenoming is a public health problem that was recently included in the list of neglected tropical diseases of the World Health Organization. In the search of new therapies for the treatment of local tissue damage induced by snake venom metalloproteinases (SVMPs), we tested the inhibitory activity of peptidomimetic compounds designed as inhibitors of matrix metalloproteinases on the activities of the SVMP Batx-I, from Bothrops atrox venom. The evaluated compounds show great potential for the inhibition of Batx-I proteolytic, hemorrhagic and edema-forming activities, especially the compound CP471474, a peptidomimetic including a hydroxamate zinc binding group. Molecular dynamics simulations suggest that binding of this compound to the enzyme is mediated by the electrostatic interaction between the hydroxamate group and the zinc cofactor, as well as contacts, mainly hydrophobic, between the side chain of the compound and amino acids located in the substrate binding subsites S1 and S1′. These results show that CP471474 constitutes a promising compound for the development of co-adjuvants to neutralize local tissue damage induced by snake venom metalloproteinases.

Taming Rugged Free Energy Landscapes Using an Average Force

The observation of complex structural transitions in biological and abiological molecular objects within time scales amenable to molecular dynamics (MD) simulations is often hampered by significant free energy barriers associated with entangled movements. Importance-sampling algorithms, a powerful class of numerical schemes for the investigation of rare events, have been widely used to extend simulations beyond the time scale common to MD. However, probing processes spanning milliseconds through microsecond molecular simulations still constitutes in practice a daunting challenge because of the difficulty of taming the ruggedness of multidimensional free energy surfaces by means of naive transition coordinates. To address this limitation, in recent years we have elaborated importance-sampling methods relying on an adaptive biasing force (ABF). In this Account, we review recent developments of algorithms aimed at mapping rugged free energy landscapes that correspond to complex processes of physical, chemical, and biological relevance. Through these developments, we have broadened the spectrum of applications of the popular ABF algorithm while improving its computational efficiency, notably for multidimensional free energy calculations. One major algorithmic advance, coined meta-eABF, merges the key features of metadynamics and an extended Lagrangian variant of ABF (eABF) by simultaneously shaving the barriers and flooding the valleys of the free energy landscape, and it possesses a convergence rate up to 5-fold greater than those of other importance-sampling algorithms. Through faster convergence and enhanced ergodic properties, meta-eABF represents a significant step forward in the simulation of millisecond-time-scale events. Here we introduce extensions of the algorithm, notably its well-tempered and replica-exchange variants, which further boost the sampling efficiency while gaining in numerical stability, thus allowing quantum-mechanical/molecular-mechanical free energy calculations to be performed at a lower cost. As a paradigm to bridge microsecond simulations to millisecond events by means of free energy calculations, we have applied the ABF family of algorithms to decompose complex movements in molecular objects of biological and abiological nature. We show here how water lubricates the shuttling of an amide-based rotaxane by altering the mechanism that underlies the concerted translation and isomerization of the macrocycle. Introducing novel collective variables in a computational workflow for the rigorous determination of standard binding free energies, we predict with utmost accuracy the thermodynamics of protein-ligand reversible association. Because of their simplicity, versatility, and robust mathematical foundations, the algorithms of the ABF family represent an appealing option for the theoretical investigation of a broad range of problems relevant to physics, chemistry, and biology.

Alchemical Free-Energy Calculations by Multiple-Replica λ-Dynamics: The Conveyor Belt Thermodynamic Integration Scheme

A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different values of the alchemical coupling parameter λ. The number K is taken to be even, and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states; and as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with K/2 equispaced replicas moving forward and K/2 equispaced replicas moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2 πK^-1, and the magnitude of its variations decreases rapidly upon increasing K, at least as K^-1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas simultaneously, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes or as a correlated multiple-replica analog of the λD or λ-local elevation umbrella sampling (λ-LEUS) schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.

Zooming across the Free-Energy Landscape: Shaving Barriers, and Flooding Valleys

A robust importance-sampling algorithm for mapping free-energy surfaces over geometrical variables, coined meta-eABF, is introduced. This algorithm shaves the free-energy barriers and floods valleys by incorporating a history-dependent potential term in the extended adaptive biasing force (eABF) framework. Numerical applications on both toy models and nontrivial examples indicate that meta-eABF explores the free-energy surface significantly faster than either eABF or metadynamics (MtD) alone, without the need to stratify the reaction pathway. In some favorable cases, meta-eABF can be as much as five times faster than other importance-sampling algorithms. Many of the shortcomings inherent to eABF and MtD, like kinetic trapping in regions of configurational space already adequately sampled, the requirement of prior knowledge of the free-energy landscape to set up the simulation, are readily eliminated in meta-eABF. Meta-eABF, therefore, represents an appealing solution for a broad range of applications, especially when both eABF and MtD fail to achieve the desired result.

The Extended Generalized Adaptive Biasing Force Algorithm for Multidimensional Free-Energy Calculations

Free-energy calculations in multiple dimensions constitute a challenging problem, owing to the significant computational cost incurred to achieve ergodic sampling. The generalized adaptive biasing force (gABF) algorithm calculates n one-dimensional lists of biasing forces to approximate the n-dimensional matrix by ignoring the coupling terms ordinarily taken into account in classical ABF simulations, thereby greatly accelerating sampling in the multidimensional space. This approximation may however occasionally lead to poor, incomplete exploration of the conformational space compared to classical ABF, especially when the selected coarse variables are strongly coupled. It has been found that introducing extended potentials coupled to the coarse variables of interest can virtually eliminate this shortcoming, and, thus, improve the efficiency of gABF simulations. In the present contribution, we propose a new free-energy method, coined extended generalized ABF (egABF), combining gABF with an extended Lagrangian strategy. The results for three illustrative examples indicate that (i) egABF can explore the transition coordinate much more efficiently compared with classical ABF, eABF, and gABF, in both simple and complex cases and (ii) egABF can achieve a higher accuracy than gABF, with a root mean-squared deviation between egABF and eABF free-energy profiles on the order of k_BT. Furthermore, the new egABF algorithm outruns the previous ABF-based algorithms in high-dimensional free-energy calculations and, hence, represents a powerful importance-sampling alternative for the investigation of complex chemical and biological processes.

Leveraging the Information from Markov State Models To Improve the Convergence of Umbrella Sampling Simulations

Umbrella sampling (US) simulation is a highly effective method for sampling the conformations of a complex system within a small subspace of predefined coordinates. In a typical US stratification strategy, biasing "window" potentials spanning the subspace of interest are introduced to narrow down the range of accessible conformations and accelerate the sampling. The speed of convergence in each biased window simulation may, however, differ. For example, windows that coincide with a large energetic barrier along a coordinate that is orthogonal to the predefined subspace are often plagued by slow relaxation timescales. Here, we design a method that can quantitatively detect this type of issue and gain further insight into the origin of the slow relaxation timescale. Once the problematic windows affected by slow convergence are identified, additional simulations limited to only these windows can be carried out, thereby reducing the overall computational effort. Several possible approaches aimed at performing US simulations adaptively are discussed, and their respective performance is illustrated using a simple model system. Last, simulations of an atomic deca-alanine system are used to demonstrate the efficacy of analyzing US simulation trajectories using the proposed method.

Extended Adaptive Biasing Force Algorithm. An On-the-Fly Implementation for Accurate Free-Energy Calculations

Proper use of the adaptive biasing force (ABF) algorithm in free-energy calculations needs certain prerequisites to be met, namely, that the Jacobian for the metric transformation and its first derivative be available and the coarse variables be independent and fully decoupled from any holonomic constraint or geometric restraint, thereby limiting singularly the field of application of the approach. The extended ABF (eABF) algorithm circumvents these intrinsic limitations by applying the time-dependent bias onto a fictitious particle coupled to the coarse variable of interest by means of a stiff spring. However, with the current implementation of eABF in the popular molecular dynamics engine NAMD, a trajectory-based post-treatment is necessary to derive the underlying free-energy change. Usually, such a posthoc analysis leads to a decrease in the reliability of the free-energy estimates due to the inevitable loss of information, as well as to a drop in efficiency, which stems from substantial read-write accesses to file systems. We have developed a user-friendly, on-the-fly code for performing eABF simulations within NAMD. In the present contribution, this code is probed in eight illustrative examples. The performance of the algorithm is compared with traditional ABF, on the one hand, and the original eABF implementation combined with a posthoc analysis, on the other hand. Our results indicate that the on-the-fly eABF algorithm (i) supplies the correct free-energy landscape in those critical cases where the coarse variables at play are coupled to either each other or to geometric restraints or holonomic constraints, (ii) greatly improves the reliability of the free-energy change, compared to the outcome of a posthoc analysis, and (iii) represents a negligible additional computational effort compared to regular ABF. Moreover, in the proposed implementation, guidelines for choosing two parameters of the eABF algorithm, namely the stiffness of the spring and the mass of the fictitious particles, are proposed. The present on-the-fly eABF implementation can be viewed as the second generation of the ABF algorithm, expected to be widely utilized in the theoretical investigation of recognition and association phenomena relevant to physics, chemistry, and biology.

Decrypting protein insertion through the translocon with free-energy calculations

Protein insertion into a membrane is a complex process involving numerous players. The most prominent of these players is the Sec translocon complex, a conserved protein-conducting channel present in the cytoplasmic membrane of bacteria and the membrane of the endoplasmic reticulum in eukaryotes. The last decade has seen tremendous leaps forward in our understanding of how insertion is managed by the translocon and its partners, coming from atomic-detailed structures, innovative experiments, and well-designed simulations. In this review, we discuss how experiments and simulations, hand-in-hand, teased out the secrets of the translocon-facilitated membrane insertion process. In particular, we focus on the role of free-energy calculations in elucidating membrane insertion. Amazingly, despite all its apparent complexity, protein insertion into membranes is primarily driven by simple thermodynamic and kinetic principles. This article is part of a Special Issue entitled: Membrane proteins edited by J.C. Gumbart and Sergei Noskov.

Why do the structural properties of complexes formed by glucans and carbon nanotubes differ so much?

Effect of glycosidic bond linkage on the structural properties of complexes formed by glucans and carbon nanotubes.

Multiple-Replica Strategies for Free-Energy Calculations in NAMD: Multiple-Walker Adaptive Biasing Force and Walker Selection Rules

From the most powerful supercomputers to multicore desktops and laptops, parallel computing architectures have been in the mainstream for some time. However, numerical schemes for calculating free energies in molecular systems that directly leverage this hardware paradigm, usually taking the form of multiple-replica strategies, are just now on the cusp of becoming standard practice. Here, we present a modification of the popular molecular dynamics program NAMD that is envisioned to facilitate the use of powerful multiple-replica strategies to improve ergodic sampling for a specific class of free-energy methods known as adaptive biasing force. We describe the software implementation in a so-called multiple-walker context, alongside the interface that makes the proposed approach accessible to the end users. We further evaluate the performance of the adaptive biasing force multiple-walker strategy for a model system, namely, the reversible folding of a short peptide, and show, in particular, in regions of the transition coordinate where convergence of the free-energy calculation is encumbered by hidden barriers, that the multiple-walker strategy can yield far more reliable results in appreciably less real time on parallel architectures, relative to standard, single-replica calculations.

Enhanced sampling techniques in molecular dynamics simulations of biological systems

BACKGROUND

Molecular dynamics has emerged as an important research methodology covering systems to the level of millions of atoms. However, insufficient sampling often limits its application. The limitation is due to rough energy landscapes, with many local minima separated by high-energy barriers, which govern the biomolecular motion.

SCOPE OF REVIEW

In the past few decades methods have been developed that address the sampling problem, such as replica-exchange molecular dynamics, metadynamics and simulated annealing. Here we present an overview over theses sampling methods in an attempt to shed light on which should be selected depending on the type of system property studied.

MAJOR CONCLUSIONS

Enhanced sampling methods have been employed for a broad range of biological systems and the choice of a suitable method is connected to biological and physical characteristics of the system, in particular system size. While metadynamics and replica-exchange molecular dynamics are the most adopted sampling methods to study biomolecular dynamics, simulated annealing is well suited to characterize very flexible systems. The use of annealing methods for a long time was restricted to simulation of small proteins; however, a variant of the method, generalized simulated annealing, can be employed at a relatively low computational cost to large macromolecular complexes.

GENERAL SIGNIFICANCE

Molecular dynamics trajectories frequently do not reach all relevant conformational substates, for example those connected with biological function, a problem that can be addressed by employing enhanced sampling algorithms. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

Thermodynamics of Deca-alanine Folding in Water

The determination of the folding dynamics of polypeptides and proteins is critical in characterizing their functions in biological systems. Numerous computational models and methods have been developed for studying structure formation at the atomic level. Due to its small size and simple structure, deca-alanine is used as a model system in molecular dynamics (MD) simulations. The free energy of unfolding in vacuum has been studied extensively using the end-to-end distance of the peptide as the reaction coordinate. However, few studies have been conducted in the presence of explicit solvent. Previous results show a significant decrease in the free energy of extended conformations in water, but the α-helical state is still notably favored over the extended state. Although sufficient in vacuum, we show that end-to-end distance is incapable of capturing the full complexity of deca-alanine folding in water. Using α-helical content as a second reaction coordinate, we deduce a more descriptive free-energy landscape, one which reveals a second energy minimum in the extended conformations that is of comparable free energy to the α-helical state. Equilibrium simulations demonstrate the relative stability of the extended and α-helical states in water as well as the transition between the two states. This work reveals both the necessity and challenge of determining a proper reaction coordinate to fully characterize a given process.