Heating and flooding: A unified approach for rapid generation of free energy surfaces

A Comprehensive Study of Factors Affecting the Prediction of the pKa Shift of Asp26 in Thioredoxin Protein

The stable protonation state of ionizable amino acids in a protein can be predicted by computing the pKa shift of that residue within the protein environment. Thermodynamic Integration (TI) is an ideal molecular dynamics-based approach for predicting the pKa shift of ionizable protein residues. Here, we probe TI-based simulation protocols for their ability to accurately predict the pKa shift of Asp26 in thioredoxin. While implicit solvent models can predict the pKa shift accurately, explicit solvent models result in substantial errors. To understand the underlying reason for this surprising discrepancy, we investigate the role of various factors such as solvent models, conformational sampling, background charges, and polarization.

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

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

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

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

Temperature Accelerated Sliced Sampling to Probe Ligand Dissociation from Protein

Modeling ligand unbinding in proteins to estimate the free energy of binding and probing the mechanism presents several challenges. They primarily pertain to the entropic bottlenecks resulting from protein and solvent conformations. While exploring the unbinding processes using enhanced sampling techniques, very long simulations are required to sample all of the conformational states as the system gets trapped in local free energy minima along transverse coordinates. Here, we demonstrate that temperature accelerated sliced sampling (TASS) is an ideal approach to overcome some of the difficulties faced by conventional sampling methods in studying ligand unbinding. Using TASS, we study the unbinding of avibactam inhibitor molecules from the Class C β-lactamase (CBL) active site. Extracting CBL-avibactam unbinding free energetics, unbinding pathways, and identifying critical interactions from the TASS simulations are demonstrated.

An exploration of machine learning models for the determination of reaction coordinates associated with conformational transitions

Determining collective variables (CVs) for conformational transitions is crucial to understanding their dynamics and targeting them in enhanced sampling simulations. Often, CVs are proposed based on intuition or prior knowledge of a system. However, the problem of systematically determining a proper reaction coordinate (RC) for a specific process in terms of a set of putative CVs can be achieved using committor analysis (CA). Identifying essential degrees of freedom that govern such transitions using CA remains elusive because of the high dimensionality of the conformational space. Various schemes exist to leverage the power of machine learning (ML) to extract an RC from CA. Here, we extend these studies and compare the ability of 17 different ML schemes to identify accurate RCs associated with conformational transitions. We tested these methods on an alanine dipeptide in vacuum and on a sarcosine dipeptoid in an implicit solvent. Our comparison revealed that the light gradient boosting machine method outperforms other methods. In order to extract key features from the models, we employed Shapley Additive exPlanations analysis and compared its interpretation with the "feature importance" approach. For the alanine dipeptide, our methodology identifies ϕ and θ dihedrals as essential degrees of freedom in the C7ax to C7eq transition. For the sarcosine dipeptoid system, the dihedrals ψ and ω are the most important for the cisαD to transαD transition. We further argue that analysis of the full dynamical pathway, and not just endpoint states, is essential for identifying key degrees of freedom governing transitions.

Exploration of high dimensional free energy landscapes by a combination of temperature-accelerated sliced sampling and parallel biasing

Temperature-accelerated sliced sampling (TASS) is an enhanced sampling method for achieving accelerated and controlled exploration of high-dimensional free energy landscapes in molecular dynamics simulations. With the aid of umbrella bias potentials, the TASS method realizes a controlled exploration and divide-and-conquer strategy for computing high-dimensional free energy surfaces. In TASS, diffusion of the system in the collective variable (CV) space is enhanced with the help of metadynamics bias and elevated-temperature of the auxiliary degrees of freedom (DOF) that are coupled to the CVs. Usually, a low-dimensional metadynamics bias is applied in TASS. In order to further improve the performance of TASS, we propose here to use a highdimensional metadynamics bias, in the same form as in a parallel bias metadynamics scheme. Here, a modified reweighting scheme, in combination with artificial neural network is used for computing unbiased probability distribution of CVs and projections of high-dimensional free energy surfaces. We first validate the accuracy and efficiency of our method in computing the four-dimensional free energy landscape for alanine tripeptide in vacuo. Subsequently, we employ the approach to calculate the eight-dimensional free energy landscape of alanine pentapeptide in vacuo. Finally, the method is applied to a more realistic problem wherein we compute the broad four-dimensional free energy surface corresponding to the deacylation of a drug molecule which is covalently complexed with a β-lactamase enzyme. We demonstrate that using parallel bias in TASS improves the efficiency of exploration of high-dimensional free energy landscapes.

Local Kernel Regression and Neural Network Approaches to the Conformational Landscapes of Oligopeptides

The application of machine learning to theoretical chemistry has made it possible to combine the accuracy of quantum chemical energetics with the thorough sampling of finite-temperature fluctuations. To reach this goal, a diverse set of methods has been proposed, ranging from simple linear models to kernel regression and highly nonlinear neural networks. Here we apply two widely different approaches to the same, challenging problem: the sampling of the conformational landscape of polypeptides at finite temperature. We develop a local kernel regression (LKR) coupled with a supervised sparsity method and compare it with a more established approach based on Behler-Parrinello type neural networks. In the context of the LKR, we discuss how the supervised selection of the reference pool of environments is crucial to achieve accurate potential energy surfaces at a competitive computational cost and leverage the locality of the model to infer which chemical environments are poorly described by the DFTB baseline. We then discuss the relative merits of the two frameworks and perform Hamiltonian-reservoir replica-exchange Monte Carlo sampling and metadynamics simulations, respectively, to demonstrate that both frameworks can achieve converged and transferable sampling of the conformational landscape of complex and flexible biomolecules with comparable accuracy and computational cost.

Boosting the conformational sampling by combining replica exchange with solute tempering and well-sliced metadynamics

Methods that combine collective variable (CV) based enhanced sampling and global tempering approaches are used in speeding-up the conformational sampling and free energy calculation of large and soft systems with a plethora of energy minima. In this paper, a new method of this kind is proposed in which the well-sliced metadynamics approach (WSMTD) is united with replica exchange with solute tempering (REST2) method. WSMTD employs a divide-and-conquer strategy wherein high-dimensional slices of a free energy surface are independently sampled and combined. The method enables one to accomplish a controlled exploration of the CV-space with a restraining bias as in umbrella sampling, and enhance-sampling of one or more orthogonal CVs using a metadynamics like bias. The new hybrid method proposed here enables boosting the sampling of more slow degrees of freedom in WSMTD simulations, without the need to specify associated CVs, through a replica exchange scheme within the framework of REST2. The high-dimensional slices of the probability distributions of CVs computed from the united WSMTD and REST2 simulations are subsequently combined using the weighted histogram analysis method to obtain the free energy surface. We show that the new method proposed here is accurate, improves the conformational sampling, and achieves quick convergence in free energy estimates. We demonstrate this by computing the conformational free energy landscapes of solvated alanine tripeptide and Trp-cage mini protein in explicit water.

Collective variable-based enhanced sampling and machine learning

ABSTRACT

Collective variable-based enhanced sampling methods have been widely used to study thermodynamic properties of complex systems. Efficiency and accuracy of these enhanced sampling methods are affected by two factors: constructing appropriate collective variables for enhanced sampling and generating accurate free energy surfaces. Recently, many machine learning techniques have been developed to improve the quality of collective variables and the accuracy of free energy surfaces. Although machine learning has achieved great successes in improving enhanced sampling methods, there are still many challenges and open questions. In this perspective, we shall review recent developments on integrating machine learning techniques and collective variable-based enhanced sampling approaches. We also discuss challenges and future research directions including generating kinetic information, exploring high-dimensional free energy surfaces, and efficiently sampling all-atom configurations.

GRAPHIC ABSTRACT

Multiscale modeling of genome organization with maximum entropy optimization

Three-dimensional (3D) organization of the human genome plays an essential role in all DNA-templated processes, including gene transcription, gene regulation, and DNA replication. Computational modeling can be an effective way of building high-resolution genome structures and improving our understanding of these molecular processes. However, it faces significant challenges as the human genome consists of over 6 × 109 base pairs, a system size that exceeds the capacity of traditional modeling approaches. In this perspective, we review the progress that has been made in modeling the human genome. Coarse-grained models parameterized to reproduce experimental data via the maximum entropy optimization algorithm serve as effective means to study genome organization at various length scales. They have provided insight into the principles of whole-genome organization and enabled de novo predictions of chromosome structures from epigenetic modifications. Applications of these models at a near-atomistic resolution further revealed physicochemical interactions that drive the phase separation of disordered proteins and dictate chromatin stability in situ. We conclude with an outlook on the opportunities and challenges in studying chromosome dynamics.

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

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

Side-Chain Polarity Modulates the Intrinsic Conformational Landscape of Model Dipeptides

The intrinsic conformational preferences of small peptides may provide additional insight into the thermodynamics and kinetics of protein folding. In this study, we explore the underlying energy landscapes of two model peptides, namely, Ac-Ala-NH₂ and Ac-Ser-NH₂, using geometry-optimization-based tools developed within the context of energy landscape theory. We analyze not only how side-chain polarity influences the structural preferences of the dipeptides, but also other emergent properties of the landscape, including heat capacity profiles, and kinetics of conformational rearrangements. The contrasting topographies of the free energy landscape agree with recent results from Fourier transform microwave spectroscopy experiments, where Ac-Ala-NH₂ was found to exist as a mixture of two conformers, while Ac-Ser-NH₂ remained structurally locked, despite exhibiting an apparently rich conformational landscape.

Molecular dynamics simulations of the effects of lipid oxidation on the permeability of cell membranes

The formation of transient pores in their membranes is a well-known mechanism of permeabilization of cells exposed to high-intensity electric pulses. However, the formation of such pores is not able to explain all aspects of the so-called electroporation phenomenon. In particular, the reasons for sustained permeability of cell membranes, persisting long after the pulses' application, remain elusive. The complete resealing of cell membranes takes indeed orders of magnitude longer than the time for electropore closure as reported from molecular dynamics (MD) investigations. Lipid peroxidation has been suggested as a possible mechanism to explain the sustainable permeability of cell membranes. However, theoretical investigations of membrane lesions containing excess amounts of hydroperoxides have shown that the conductivities of such lesions were not high enough to account for the experimental measurements. Here, expanding on these studies, we investigate quantitatively the permeability of cell membrane lesions that underwent secondary oxidation. MD simulations and free energy calculations of lipid bilayers show that such lesions provide a better model of post-pulse permeable and conductive electropermeabilized cells. These results are further discussed in the context of sonoporation and ferroptosis, respectively a procedure and a phenomenon, among others, in which, alike electroporation, substantial lipid oxidation might be triggered.

Stability and folding pathways of tetra-nucleosome from six-dimensional free energy surface

The three-dimensional organization of chromatin is expected to play critical roles in regulating genome functions. High-resolution characterization of its structure and dynamics could improve our understanding of gene regulation mechanisms but has remained challenging. Using a near-atomistic model that preserves the chemical specificity of protein-DNA interactions at residue and base-pair resolution, we studied the stability and folding pathways of a tetra-nucleosome. Dynamical simulations performed with an advanced sampling technique uncovered multiple pathways that connect open chromatin configurations with the zigzag crystal structure. Intermediate states along the simulated folding pathways resemble chromatin configurations reported from in situ experiments. We further determined a six-dimensional free energy surface as a function of the inter-nucleosome distances via a deep learning approach. The zigzag structure can indeed be seen as the global minimum of the surface. However, it is not favored by a significant amount relative to the partially unfolded, in situ configurations. Chemical perturbations such as histone H4 tail acetylation and thermal fluctuations can further tilt the energetic balance to stabilize intermediate states. Our study provides insight into the connection between various reported chromatin configurations and has implications on the in situ relevance of the 30 nm fiber.

The three-dimensional organization of chromatin plays critical roles in regulating genome function. Here the authors apply a near atomistic model to study the structure and dynamics of the chromatin folding unit - the tetra-nucleosome - to provide insight into how chromatin folds.

Refined metadynamics through canonical sampling using time-invariant bias potential: A study of polyalcohol dehydration in hot acidic solutions

We propose a canonical sampling method to refine metadynamics simulations a posteriori, where the hills obtained from metadynamics are used as a time-invariant bias potential. In this way, the statistical error in the computed reaction barriers is reduced by an efficient sampling of the collective variable space at the free energy level of interest. This simple approach could be useful particularly when two or more free energy barriers are to be compared among chemical reactions in different or competing conditions. The method was then applied to study the acid dependence of polyalcohol dehydration reactions in high-temperature aqueous solutions. It was found that the reaction proceeds consistently via an S_N 2 mechanism, whereby the free energy of protonation of the hydroxyl group created as an intermediate is affected significantly by the acidic species. Although demonstration is shown for a specific problem, the computational method suggested herein could be generally used for simulations of complex reactions in the condensed phase.

New insights on the role of ROS in the mechanisms of sonoporation-mediated gene delivery

Reactive oxygen species (ROS) are hypothesized to play a role in the sonoporation mechanisms. Nevertheless, the acoustical phenomenon behind the ROS production as well as the exact mechanisms of ROS action involved in the increased cell membrane permeability are still not fully understood. Therefore, we investigated the key processes occurring at the molecular level in and around microbubbles subjected to ultrasound using computational chemistry methods. To confirm the molecular simulation predictions, we measured the ROS production by exposing SonoVue® microbubbles (MBs) to ultrasound using biological assays. To investigate the role of ROS in cell membrane permeabilization, cells were subjected to ultrasound in presence of MBs and plasmid encoding reporter gene, and the transfection level was assessed using flow cytometry. The molecular simulations showed that under sonoporation conditions, ROS can form inside the MBs. These radicals could easily diffuse through the MB shell toward the surrounding aqueous phase and participate in the permeabilization of nearby cell membranes. Experimental data confirmed that MBs favor spontaneous formation of a host of free radicals where HO was the main ROS species after US exposure. The presence of ROS scavengers/inhibitors during the sonoporation process decreased both the production of ROS and the subsequent transfection level without significant loss of cell viability. In conclusion, the exposure of MBs to ultrasound might be the origin of chemical effects, which play a role in the cell membrane permeabilization and in the in vitro gene delivery when generated in its proximity.

Analytical Method Using a Scaled Hypersphere Search for High-Dimensional Metadynamics Simulations

Metadynamics (MTD) is one of the most effective methods for calculating the free energy surface and finding rare events. Nevertheless, numerous studies using MTD have been carried out using 3D or lower dimensional collective variables (CVs), as higher dimensional CVs require costly computational resources and the obtained results are too complex to understand the important events. The latter issue can be conveniently solved by utilizing the free energy reaction network (FERN), which is a graph structure consisting of edges of minimum free energy paths (MFEPs), nodes of equation (EQ) points, and transition state (TS) points. In the present article, a new method for exploring FERNs on high-dimensional CVs using MTD and the scaled hypersphere search (SHS) method is described. A test calculation based on the MTD-SHS simulation of met-enkephalin in explicit water with 7 CVs was conducted. As a result, 889 EQ points and 1805 TS points were found. The MTD-SHS approach can find MFEPs exhaustively; therefore, the FERNs can be estimated without any a priori knowledge of the EQ and TS points.

Comparison of the Performance of Machine Learning Models in Representing High-Dimensional Free Energy Surfaces and Generating Observables

Free energy surfaces of chemical and physical systems are often generated using a popular class of enhanced sampling methods that target a set of collective variables (CVs) chosen to distinguish the characteristic features of these surfaces. While some of these approaches are typically limited to low (∼1-3)-dimensional CV subspaces, methods such as driven adiabatic free-energy dynamics/temperature-accelerated molecular dynamics have been shown to be capable of generating free energy surfaces of quite high dimension by sampling the associated marginal probability distribution via full sweeps over the CV landscape. These approaches repeatedly visit conformational basins, producing a scattering of points within the basins on each visit. Consequently, they are particularly amenable to synergistic combination with regression machine learning methods for filling in the surfaces between the sampled points and for providing a compact and continuous (or semicontinuous) representation of the surfaces that can be easily stored and used for further computation of observable properties. Given the central role of machine learning techniques in this combined approach, it is timely to provide a detailed comparison of the performance of different machine learning strategies and models, including neural networks, kernel ridge regression, support vector machines, and weighted neighbor schemes, for their ability to learn these high-dimensional surfaces as a function of the amount of sampled training data and, once trained, to subsequently generate accurate ensemble averages corresponding to observable properties of the systems. In this article, we perform such a comparison on a set of oligopeptides, in both gas and aqueous phases, corresponding to CV spaces of 2-10 dimensions and assess their ability to provide a global representation of the free energy surfaces and to generate accurate ensemble averages.

Neural-Network-Based Path Collective Variables for Enhanced Sampling of Phase Transformations

The investigation of the microscopic processes underlying structural phase transformations in solids is extremely challenging for both simulation and experiment. Atomistic simulations of solid-solid phase transitions require extensive sampling of the corresponding high-dimensional and often rugged energy landscape. Here, we propose a rigorous construction of a 1D path collective variable that is used in combination with enhanced sampling techniques for efficient exploration of the transformation mechanisms. The path collective variable is defined in a space spanned by global classifiers that are derived from local structural units. A reliable identification of the local structural environments is achieved by employing a neural-network-based classification scheme. The proposed path collective variable is generally applicable and enables the investigation of both transformation mechanisms and kinetics.

Mapping saddles and minima on free energy surfaces using multiple climbing strings

Locating saddle points on free energy surfaces is key in characterizing multistate transition events in complicated molecular-scale systems. Because these saddle points represent transition states, determining minimum free energy pathways to these saddles and measuring their free energies relative to their connected minima are further necessary, for instance, to estimate transition rates. In this work, we propose a new multistring version of the climbing string method in collective variables to locate all saddles and corresponding pathways on free energy surfaces. The method uses dynamic strings to locate saddles and static strings to keep a history of prior strings converged to saddles. Interaction of the dynamic strings with the static strings is used to avoid the convergence to already-identified saddles. Additionally, because the strings approximate curves in collective-variable space, and we can measure free energy along each curve, identification of any saddle's two connected minima is guaranteed. We demonstrate this method to map the network of stationary points in the 2D and 4D free energy surfaces of alanine dipeptide and alanine tripeptide, respectively.

Exploration of Chemical Compound, Conformer, and Reaction Space with Meta-Dynamics Simulations Based on Tight-Binding Quantum Chemical Calculations

The semiempirical tight-binding based quantum chemistry method GFN2-xTB is used in the framework of meta-dynamics (MTD) to globally explore chemical compound, conformer, and reaction space. The biasing potential given as a sum of Gaussian functions is expressed with the root-mean-square-deviation (RMSD) in Cartesian space as a metric for the collective variables. This choice makes the approach robust and generally applicable to three common problems (i.e., conformer search, chemical reaction space exploration in a virtual nanoreactor, and for guessing reaction paths). Because of the inherent locality of the atomic RMSD, functional group or fragment selective treatments are possible facilitating the investigation of catalytic processes where, for example, only the substrate is thermally activated. Due to the approximate character of the GFN2-xTB method, the resulting structure ensembles require further refinement with more sophisticated, for example, density functional or wave function theory methods. However, the approach is extremely efficient running routinely on common laptop computers in minutes to hours of computation time even for realistically sized molecules with a few hundred atoms. Furthermore, the underlying potential energy surface for molecules containing almost all elements ( Z = 1-86) is globally consistent including the covalent dissociation process and electronically complicated situations in, for example, transition metal systems. As examples, thermal decomposition, ethyne oligomerization, the oxidation of hydrocarbons (by oxygen and a P450 enzyme model), a Miller-Urey model system, a thermally forbidden dimerization, and a multistep intramolecular cyclization reaction are shown. For typical conformational search problems of organic drug molecules, the new MTD(RMSD) algorithm yields lower energy structures and more complete conformer ensembles at reduced computational effort compared with its already well performing predecessor.

Finding Free-Energy Landmarks of Chemical Reactions

We propose a novel approach to search for free-energy landmarks, i.e., minima and the saddle points, of chemical reactions in an automated manner using a combination of steepest descent and gentlest ascent methods. A numerical approach is suggested to improve the sampling efficiency of the second derivatives of the free-energy surface, which is required in the gentlest ascent method. This technique opens a way to identify free-energy landmarks of bond-breaking/creating processes in which the underlying potential energy surface is described using on-the-fly electronic structure calculations. As demonstrations of the approach, we present applications to the ring-opening of cis-1,2-dimethylbenzocyclobutene using the semiempirical PM7 method, focusing on the temperature dependence of the paths and barrier of the reaction, and we study an SN2 reaction in aqueous solution using a semiempirical QM/MM approach combining PM7 with the TIP3P water model.

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.

Molecular dynamics based enhanced sampling of collective variables with very large time steps

Enhanced sampling techniques that target a set of collective variables and that use molecular dynamics as the driving engine have seen widespread application in the computational molecular sciences as a means to explore the free-energy landscapes of complex systems. The use of molecular dynamics as the fundamental driver of the sampling requires the introduction of a time step whose magnitude is limited by the fastest motions in a system. While standard multiple time-stepping methods allow larger time steps to be employed for the slower and computationally more expensive forces, the maximum achievable increase in time step is limited by resonance phenomena, which inextricably couple fast and slow motions. Recently, we introduced deterministic and stochastic resonance-free multiple time step algorithms for molecular dynamics that solve this resonance problem and allow ten- to twenty-fold gains in the large time step compared to standard multiple time step algorithms [P. Minary et al., Phys. Rev. Lett. 93, 150201 (2004); B. Leimkuhler et al., Mol. Phys. 111, 3579-3594 (2013)]. These methods are based on the imposition of isokinetic constraints that couple the physical system to Nosé-Hoover chains or Nosé-Hoover Langevin schemes. In this paper, we show how to adapt these methods for collective variable-based enhanced sampling techniques, specifically adiabatic free-energy dynamics/temperature-accelerated molecular dynamics, unified free-energy dynamics, and by extension, metadynamics, thus allowing simulations employing these methods to employ similarly very large time steps. The combination of resonance-free multiple time step integrators with free-energy-based enhanced sampling significantly improves the efficiency of conformational exploration.

Effects of hydration on the protonation state of a lysine analog crossing a phospholipid bilayer – insights from molecular dynamics and free-energy calculations

Protonation states of amino acids crossing lipid bilayers from multidimensional free energy surfaces.

Stochastic Neural Network Approach for Learning High-Dimensional Free Energy Surfaces

The generation of free energy landscapes corresponding to conformational equilibria in complex molecular systems remains a significant computational challenge. Adding to this challenge is the need to represent, store, and manipulate the often high-dimensional surfaces that result from rare-event sampling approaches employed to compute them. In this Letter, we propose the use of artificial neural networks as a solution to these issues. Using specific examples, we discuss network training using enhanced-sampling methods and the use of the networks in the calculation of ensemble averages.

Optimization and Automation of the Construction of Smooth Free Energy Profiles

An adaptive procedure is introduced to construct smooth analytical profiles of the free energy along a reaction coordinate using sampled data from multiple biased simulations. The procedure is based upon identifying problematic regions encountered in maximum likelihood estimators of the profile where there are statistically relevant discrepancies between the empirical and parametrized cumulative distribution functions and preferentially improving the construction of the parametric profile in these regions. The method is designed to produce continuous and smooth analytical fits that satisfy statistical goodness-of-fit tests with a minimum number of parameters. The accuracy of the profile obtained from the adaptive construction is compared by numerical computation to that of smooth interpolations based on an optimally chosen weighted histogram method for a solvated ion pair system and for an activated process for which the analytical form of the potential of mean force is available. In the model where the exact profile is known, the adaptive procedure is shown to reduce the integrated error relative to the optimal histogram construction by a factor of 3 or more in the typical case where the sampling is not extensive. It is demonstrated that the adaptive procedure can be used to produce statistically accurate smooth analytical representations of the free energy profile that can be evaluated with little computational effort and require little user input.

Reliable and efficient reaction path and transition state finding for surface reactions with the growing string method

The computational challenge of fast and reliable transition state and reaction path optimization requires new methodological strategies to maintain low cost, high accuracy, and systematic searching capabilities. The growing string method using internal coordinates has proven to be highly effective for the study of molecular, gas phase reactions, but difficulties in choosing a suitable coordinate system for periodic systems has prevented its use for surface chemistry. New developments are therefore needed, and presented herein, to handle surface reactions which include atoms with large coordination numbers that cannot be treated using standard internal coordinates. The double-ended and single-ended growing string methods are implemented using a hybrid coordinate system, then benchmarked for a test set of 43 elementary reactions occurring on surfaces. These results show that the growing string method is at least 45% faster than the widely used climbing image-nudged elastic band method, which also fails to converge in several of the test cases. Additionally, the surface growing string method has a unique single-ended search method which can move outward from an initial structure to find the intermediates, transition states, and reaction paths simultaneously. This powerful explorative feature of single ended-growing string method is demonstrated to uncover, for the first time, the mechanism for atomic layer deposition of TiN on Cu(111) surface. This reaction is found to proceed through multiple hydrogen-transfer and ligand-exchange events, while formation of H-bonds stabilizes intermediates of the reaction. Purging gaseous products out of the reaction environment is the driving force for these reactions. © 2017 Wiley Periodicals, Inc.

Effective Riemannian Diffusion Model for Conformational Dynamics of Biomolecular Systems

We present a Riemannian formalism for effective diffusion of biomolecules in collective variable spaces that provides a robust framework for conformational free energy calculation methods. Unlike their Euclidean counterparts, the Riemannian potential of mean force (PMF) and minimum free energy path (MFEP) are invariant under coordinate transformations. The presented formalism can be readily employed to modify the collective variable based enhanced sampling techniques, such as umbrella sampling (US) commonly used in biomolecular simulations, to take into account the role of intrinsic geometry of collective variable space. Although our model is mathematically equivalent to a Euclidean diffusion with a position-dependent diffusion tensor, the Riemannian formulation provides a more convenient framework for free energy calculation methods and path-finding algorithms aimed at characterizing the effective conformational dynamics of biomolecules. A simple three-dimensional toy model and a pentapeptide (met-enkephalin) simulated in an explicit solvent environment are used to illustrate the workings of the formalism and its implementation.

Molecular dynamics simulations of cholesterol-rich membranes using a coarse-grained force field for cyclic alkanes

The architecture of a biological membrane hinges upon the fundamental fact that its properties are determined by more than the sum of its individual components. Studies on model membranes have shown the need to characterize in molecular detail how properties such as thickness, fluidity, and macroscopic bending rigidity are regulated by the interactions between individual molecules in a non-trivial fashion. Simulation-based approaches are invaluable to this purpose but are typically limited to short sampling times and model systems that are often smaller than the required properties. To alleviate both limitations, the use of coarse-grained (CG) models is nowadays an established computational strategy. We here present a new CG force field for cholesterol, which was developed by using measured properties of small molecules, and can be used in combination with our previously developed force field for phospholipids. The new model performs with precision comparable to atomistic force fields in predicting the properties of cholesterol-rich phospholipid bilayers, including area per lipid, bilayer thickness, tail order parameter, increase in bending rigidity, and propensity to form liquid-ordered domains in ternary mixtures. We suggest the use of this model to quantify the impact of cholesterol on macroscopic properties and on microscopic phenomena involving localization and trafficking of lipids and proteins on cellular membranes.

A Stochastic, Resonance-Free Multiple Time-Step Algorithm for Polarizable Models That Permits Very Large Time Steps

Molecular dynamics remains one of the most widely used computational tools in the theoretical molecular sciences to sample an equilibrium ensemble distribution and/or to study the dynamical properties of a system. The efficiency of a molecular dynamics calculation is limited by the size of the time step that can be employed, which is dictated by the highest frequencies in the system. However, many properties of interest are connected to low-frequency, long time-scale phenomena, requiring many small time steps to capture. This ubiquitous problem can be ameliorated by employing multiple time-step algorithms, which assign different time steps to forces acting on different time scales. In such a scheme, fast forces are evaluated more frequently than slow forces, and as the former are often computationally much cheaper to evaluate, the savings can be significant. Standard multiple time-step approaches are limited, however, by resonance phenomena, wherein motion on the fastest time scales limits the step sizes that can be chosen for the slower time scales. In atomistic models of biomolecular systems, for example, the largest time step is typically limited to around 5 fs. Previously, we introduced an isokinetic extended phase-space algorithm (Minary et al. Phys. Rev. Lett. 2004, 93, 150201) and its stochastic analog (Leimkuhler et al. Mol. Phys. 2013, 111, 3579) that eliminate resonance phenomena through a set of kinetic energy constraints. In simulations of a fixed-charge flexible model of liquid water, for example, the time step that could be assigned to the slow forces approached 100 fs. In this paper, we develop a stochastic isokinetic algorithm for multiple time-step molecular dynamics calculations using a polarizable model based on fluctuating dipoles. The scheme developed here employs two sets of induced dipole moments, specifically, those associated with short-range interactions and those associated with a full set of interactions. The scheme is demonstrated on the polarizable AMOEBA water model. As was seen with fixed-charge models, we are able to obtain large time steps exceeding 100 fs, allowing calculations to be performed 10 to 20 times faster than standard thermostated molecular dynamics.

Locating landmarks on high-dimensional free energy surfaces

Coarse graining of complex systems possessing many degrees of freedom can often be a useful approach for analyzing and understanding key features of these systems in terms of just a few variables. The relevant energy landscape in a coarse-grained description is the free energy surface as a function of the coarse-grained variables, which, despite the dimensional reduction, can still be an object of high dimension. Consequently, navigating and exploring this high-dimensional free energy surface is a nontrivial task. In this paper, we use techniques from multiscale modeling, stochastic optimization, and machine learning to devise a strategy for locating minima and saddle points (termed "landmarks") on a high-dimensional free energy surface "on the fly" and without requiring prior knowledge of or an explicit form for the surface. In addition, we propose a compact graph representation of the landmarks and connections between them, and we show that the graph nodes can be subsequently analyzed and clustered based on key attributes that elucidate important properties of the system. Finally, we show that knowledge of landmark locations allows for the efficient determination of their relative free energies via enhanced sampling techniques.

Practical Aspects of Free-Energy Calculations: A Review

Free-energy calculations in the framework of classical molecular dynamics simulations are nowadays used in a wide range of research areas including solvation thermodynamics, molecular recognition, and protein folding. The basic components of a free-energy calculation, that is, a suitable model Hamiltonian, a sampling protocol, and an estimator for the free energy, are independent of the specific application. However, the attention that one has to pay to these components depends considerably on the specific application. Here, we review six different areas of application and discuss the relative importance of the three main components to provide the reader with an organigram and to make nonexperts aware of the many pitfalls present in free energy calculations.