Estimation Of Free Energy Differences Research Articles

How to Sample Dozens of Substitutions per Site with λ Dynamics.

Alchemical free energy methods are useful in computer-aided drug design and computational protein design because they provide rigorous statistical mechanics-based estimates of free energy differences from molecular dynamics simulations. λ dynamics is a free energy method with the ability to characterize combinatorial chemical spaces spanning thousands of related systems within a single simulation, which gives it a distinct advantage over other alchemical free energy methods that are mostly limited to pairwise comparisons. Recently developed methods have improved the scalability of λ dynamics to perturbations at many sites; however, the size of chemical space that can be explored at each individual site has previously been limited to fewer than ten substituents. As the number of substituents increases, the volume of alchemical space corresponding to nonphysical alchemical intermediates grows exponentially relative to the size corresponding to the physical states of interest. Beyond nine substituents, λ dynamics simulations become lost in an alchemical morass of intermediate states. In this work, we introduce new biasing potentials that circumvent excessive sampling of intermediate states by favoring sampling of physical end points relative to alchemical intermediates. Additionally, we present a more scalable adaptive landscape flattening algorithm for these larger alchemical spaces. Finally, we show that this potential enables more efficient sampling in both protein and drug design test systems with up to 24 substituents per site, enabling, for the first time, simultaneous simulation of all 20 amino acids.

Force Field X: A computational microscope to study genetic variation and organic crystals using theory and experiment.

Force Field X (FFX) is an open-source software package for atomic resolution modeling of genetic variants and organic crystals that leverages advanced potential energy functions and experimental data. FFX currently consists of nine modular packages with novel algorithms that include global optimization via a many-body expansion, acid-base chemistry using polarizable constant-pH molecular dynamics, estimation of free energy differences, generalized Kirkwood implicit solvent models, and many more. Applications of FFX focus on the use and development of a crystal structure prediction pipeline, biomolecular structure refinement against experimental datasets, and estimation of the thermodynamic effects of genetic variants on both proteins and nucleic acids. The use of Parallel Java and OpenMM combines to offer shared memory, message passing, and graphics processing unit parallelization for high performance simulations. Overall, the FFX platform serves as a computational microscope to study systems ranging from organic crystals to solvated biomolecular systems.

On the calculation of free energies over Hamiltonian and order parameters via perturbation and thermodynamic integration.

In this work, complementary formulas are presented to compute free-energy differences via perturbation (FEP) methods and thermodynamic integration (TI). These formulas are derived by selecting only the most statistically significant data from the information extractable from the simulated points involved. On the one hand, commonly used FEP techniques based on overlap sampling leverage the full information contained in the overlapping macrostate probability distributions. On the other hand, conventional TI methods only use information on the first moments of those distributions, as embodied by the first derivatives of the free energy. Since the accuracy of simulation data degrades considerably for high-order moments (for FEP) or free-energy derivatives (for TI), it is proposed to consider, consistently for both methods, data up to second-order moments/derivatives. This provides a compromise between the limiting strategies embodied by common FEP and TI and leads to simple, optimized expressions to evaluate free-energy differences. The proposed formulas are validated with an analytically solvable harmonic Hamiltonian (for assessing systematic errors), an atomistic system (for computing the potential of mean force with coordinate-dependent order parameters), and a binary-component coarse-grained model (for tracing a solid-liquid phase diagram in an ensemble sampled through alchemical transformations). It is shown that the proposed FEP and TI formulas are straightforward to implement, perform similarly well, and allow robust estimation of free-energy differences even when the spacing of successive points does not guarantee them to have proper overlapping in phase space.

Automated, Accurate, and Scalable Relative Protein-Ligand Binding Free-Energy Calculations Using Lambda Dynamics.

Accurate predictions of changes to protein-ligand binding affinity in response to chemical modifications are of utility in small-molecule lead optimization. Relative free-energy perturbation (FEP) approaches are one of the most widely utilized for this goal but involve significant computational cost, thus limiting their application to small sets of compounds. Lambda dynamics, also rigorously based on the principles of statistical mechanics, provides a more efficient alternative. In this paper, we describe the development of a workflow to set up, execute, and analyze multisite lambda dynamics (MSLD) calculations run on GPUs with CHARMM implemented in BIOVIA Discovery Studio and Pipeline Pilot. The workflow establishes a framework for setting up simulation systems for exploratory screening of modifications to a lead compound, enabling the calculation of relative binding affinities of combinatorial libraries. To validate the workflow, a diverse data set of congeneric ligands for seven proteins with experimental binding affinity data is examined. A protocol to automatically tailor fit biasing potentials iteratively to flatten the free-energy landscape of any MSLD system is developed, which enhances sampling and allows for efficient estimation of free-energy differences. The protocol is first validated on a large number of ligand subsets that model diverse substituents, which shows accurate and reliable performance. The scalability of the workflow is also tested to screen more than 100 ligands modeled in a single system, which also resulted in accurate predictions. With a cumulative sampling time of 150 ns or less, the method results in average unsigned errors of under 1 kcal/mol in most cases for both small and large combinatorial libraries. For the multisite systems examined, the method is estimated to be more than an order of magnitude more efficient than contemporary FEP applications. The results thus demonstrate the utility of the presented MSLD workflow to efficiently screen combinatorial libraries and explore the chemical space around a lead compound and thus are of utility in lead optimization.

Ligand Binding Thermodynamic Cycles: Hysteresis, the Locally Weighted Histogram Analysis Method, and the Overlapping States Matrix.

Free energy perturbation (FEP) simulations have been widely applied to obtain predictions of the relative binding free energy for a series of congeneric ligands binding to the same receptor, which is an essential component for the lead optimization process in computer-aided drug discovery. In the case of several congeneric ligands forming a perturbation map involving a closed thermodynamic cycle, the summation of the estimated free energy change along each edge in the cycle using Bennett acceptance ratio (BAR) usually will deviate from zero due to systematic and random errors, which is the hysteresis of cycle closure. In this work, the advanced reweighting techniques binless weighted histogram analysis method (UWHAM) and locally weighted histogram analysis method (LWHAM) are applied to provide statistical estimators of the free energy change along each edge in order to eliminate the hysteresis effect. As an example, we analyze a closed thermodynamic cycle involving four congeneric ligands which bind to HIV-1 integrase, a promising target which has emerged for antiviral therapy. We demonstrate that, compared with FEP and BAR, more accurate and hysteresis-free estimates of free energy differences can be achieved by using UWHAM to find a single estimate of the density of states based on all of the data in the cycle. Furthermore, by comparison of LWHAM results obtained from the inclusion of different numbers of neighboring states with UWHAM estimation involving all the states, we show how to determine the optimal neighborhood size in the LWHAM analysis to balance the trade-offs between computational cost and accuracy of the free energy prediction. Even with the smallest neighborhood, LWHAM can improve the BAR free energy estimates using the same input data as BAR. We introduce an overlapping states matrix that is constructed by using the global jump formula of LWHAM and plot its heat map. The heat map provides a quantitative measure of the overlap between pairs of alchemical/thermodynamic states. We explain how to identify and improve the FEP calculations along the edges that most likely cause large systematic errors by using the heat map of the overlapping states matrix and by comparing the BAR and UWHAM estimates of the free energy change.

Molecular system identification for enzyme directed evolution and design.

The rational design of chemical catalysts requires methods for the measurement of free energy differences in the catalytic mechanism for any given catalyst Hamiltonian. The scope of experimental learning algorithms that can be applied to catalyst design would also be expanded by the availability of such methods. Methods for catalyst characterization typically either estimate apparent kinetic parameters that do not necessarily correspond to free energy differences in the catalytic mechanism or measure individual free energy differences that are not sufficient for establishing the relationship between the potential energy surface and catalytic activity. Moreover, in order to enhance the duty cycle of catalyst design, statistically efficient methods for the estimation of the complete set of free energy differences relevant to the catalytic activity based on high-throughput measurements are preferred. In this paper, we present a theoretical and algorithmic system identification framework for the optimal estimation of free energy differences in solution phase catalysts, with a focus on one- and two-substrate enzymes. This framework, which can be automated using programmable logic, prescribes a choice of feasible experimental measurements and manipulated input variables that identify the complete set of free energy differences relevant to the catalytic activity and minimize the uncertainty in these free energy estimates for each successive Hamiltonian design. The framework also employs decision-theoretic logic to determine when model reduction can be applied to improve the duty cycle of high-throughput catalyst design. Automation of the algorithm using fluidic control systems is proposed, and applications of the framework to the problem of enzyme design are discussed.

Efficient Round-Trip Time Optimization for Replica-Exchange Enveloping Distribution Sampling (RE-EDS)

Replica-exchange enveloping distribution sampling (RE-EDS) allows the efficient estimation of free-energy differences between multiple end-states from a single molecular dynamics (MD) simulation. In EDS, a reference state is sampled, which can be tuned by two types of parameters, i.e., smoothness parameters(s) and energy offsets, such that all end-states are sufficiently sampled. However, the choice of these parameters is not trivial. Replica exchange (RE) or parallel tempering is a widely applied technique to enhance sampling. By combining EDS with the RE technique, the parameter choice problem could be simplified and the challenge shifted toward an optimal distribution of the replicas in the smoothness-parameter space. The choice of a certain replica distribution can alter the sampling efficiency significantly. In this work, global round-trip time optimization (GRTO) algorithms are tested for the use in RE-EDS simulations. In addition, a local round-trip time optimization (LRTO) algorithm is proposed for systems with slowly adapting environments, where a reliable estimate for the round-trip time is challenging to obtain. The optimization algorithms were applied to RE-EDS simulations of a system of nine small-molecule inhibitors of phenylethanolamine N-methyltransferase (PNMT). The energy offsets were determined using our recently proposed parallel energy-offset (PEOE) estimation scheme. While the multistate GRTO algorithm yielded the best replica distribution for the ligands in water, the multistate LRTO algorithm was found to be the method of choice for the ligands in complex with PNMT. With this, the 36 alchemical free-energy differences between the nine ligands were calculated successfully from a single RE-EDS simulation 10 ns in length. Thus, RE-EDS presents an efficient method for the estimation of relative binding free energies.

A Monte Carlo Resampling Approach for the Calculation of Hybrid Classical and Quantum Free Energies.

Hybrid free energy methods allow estimation of free energy differences at the quantum mechanics (QM) level with high efficiency by performing sampling at the classical mechanics (MM) level. Various approaches to allow the calculation of QM corrections to classical free energies have been proposed. The single step free energy perturbation approach starts with a classically generated ensemble, a subset of structures of which are postprocessed to obtain QM energies for use with the Zwanzig equation. This gives an estimate of the free energy difference associated with the change from an MM to a QM Hamiltonian. Owing to the poor numerical properties of the Zwanzig equation, however, recent developments have produced alternative methods which aim to provide access to the properties of the true QM ensemble. Here we propose an approach based on the resampling of MM structural ensembles and application of a Monte Carlo acceptance test which in principle, can generate the exact QM ensemble or intermediate ensembles between the MM and QM states. We carry out a detailed comparison against the Zwanzig equation and recently proposed non-Boltzmann methods. As a test system we use a set of small molecule hydration free energies for which hybrid free energy calculations are performed at the semiempirical Density Functional Tight Binding level. Equivalent ensembles at this level of theory have also been generated allowing the reverse QM to MM perturbations to be performed along with a detailed analysis of the results. Additionally, a previously published nucleotide base pair data set simulated at the QM level using ab initio molecular dynamics is also considered. We provide a strong rationale for the use of the Monte Carlo Resampling and non-Boltzmann approaches by showing that configuration space overlaps can be estimated which provide useful diagnostic information regarding the accuracy of these hybrid approaches.

Symmetry for the duration of entropy-consuming intervals.

We introduce the violation fraction υ as the cumulative fraction of time that a mesoscopic system spends consuming entropy at a single trajectory in phase space. We show that the fluctuations of this quantity are described in terms of a symmetry relation reminiscent of fluctuation theorems, which involve a function Φ, which can be interpreted as an entropy associated with the fluctuations of the violation fraction. The function Φ, when evaluated for arbitrary stochastic realizations of the violation fraction, is odd upon the symmetry transformations that are relevant for the associated stochastic entropy production. This fact leads to a detailed fluctuation theorem for the probability density function of Φ. We study the steady-state limit of this symmetry in the paradigmatic case of a colloidal particle dragged by optical tweezers through an aqueous solution. Finally, we briefly discuss possible applications of our results for the estimation of free-energy differences from single-molecule experiments.

Analysis of Free Energy of Opening the Troponin C Binding Pocket for Troponin I using Microsecond Molecular Dynamics Simulations

Troponin (Tn), part of the thin filament in cardiomyocytes, plays an important role in calcium signaling events in cardiac muscle contraction. It acts as a Ca2+-dependent switch, activating and deactivating the myofilament leading to contraction and relaxation of the muscle cell. An important mechanism in the regulation of contraction is the opening up of the TnC hydrophobic patch to allow TnI to bind. We performed microsecond molecular dynamics simulations of TnC and two of its mutants (V44Q, E40A) in different states of calcium binding. The dynamics and kinetics of the observed opening transitions allow for the estimation of free energy differences for opening. We obtained values for wildtype Ca2+-bound TnC (∼8 kcal/mol), V44Q Ca2+-bound TnC (3.2 kcal/mol), E40A Ca2+-bound TnC (∼12 kcal/mol) and wildtype apo TnC (∼20 kcal/mol). These results suggest that the mutations have profound impact on the frequency of presenting the hydrophobic patch to TnI. These simulations also corroborate that cardiac wildtype TnC does not open on timescales relevant to contraction without calcium being bound. Additionally, Brownian dynamics simulations are used to investigate TnI association with TnC. Simulations of a full length troponin model elucidate the dynamical interplay between the TnC, TnI and TnT subunits.

Estimation of Free-Energy Differences from Computed Work Distributions: An Application of Jarzynski’s Equality

Equilibrium free-energy differences can be computed from nonequilibrium molecular dynamics (MD) simulations using Jarzynski's equality (Jarzynski, C. Phys. Rev. Lett.1997, 78, 2690) by combining a large set of independent trajectories (path ensemble). Here we present the multistep trajectory combination (MSTC) method to compute free-energy differences, which by combining trajectories significantly reduces the number of trajectories necessary to generate a representative path ensemble. This method generates well-sampled work distributions, even for large systems, by combining parts of a relatively small number of trajectories carried out in steps. To assess the efficiency of the MSTC method, we derived analytical expressions and used them to compute the bias and the variance of the free-energy estimates along with numerically calculated values. We show that the MSTC method significantly reduces both the bias and variance of the free-energy estimates compared to the estimates obtained using single-step trajectories. In addition, because in the MSTC method the process is divided into steps, it is feasible to compute the reverse transition. By combining the forward and reverse processes, the free-energy difference can be computed using the Crooks' fluctuation theorem (Crooks, G. E. J. Stat. Phys.1998, 90, 1481 and Crooks, G. E. Phys. Rev. E 2000, 61, 2361) or Bennett's acceptance ratio (Bennett, C. H. J. Comput. Phys. 1976, 22, 245), which further reduces the bias and variance of the estimates.

Linear‐scaling soft‐core scheme for alchemical free energy calculations

Alchemical free energy calculations involving the removal or insertion of atoms into condensed phase systems generally make use of soft-core scaling of nonbonded interactions, designed to circumvent numerical instabilities that arise from weakly interacting "hard" atoms in close proximity. Current methods model soft-core atoms by introducing a nonlinear dependence between the shape of the interaction potential and the strength of the interaction. In this article, we propose a soft-core method that avoids introducing such a nonlinear dependence, through the application of a smooth flattening of the potential energy only in a region that is energetically accessible under normal conditions. We discuss the benefits that this entails and explore a selection of applications, including enhanced methods for the estimation of free energy differences and for the automated optimization of the placement of intermediate states in multistage alchemical calculations.

Osmotic ensemble methods for predicting adsorption-induced structural transitions in nanoporous materials using molecular simulations

Osmotic framework adsorbed solution theory is a useful molecular simulation method to predict the evolution of structural transitions upon adsorption of guest molecules in flexible nanoporous solids. One challenge with previous uses of this approach has been the estimation of free energy differences between the solid phases of interest in the absence of adsorbed molecules. Here we demonstrate that these free energy differences can be calculated without reference to experimental data via the vibrational density of states of each phase, a quantity that can be obtained from molecular dynamics simulations. We show the applicability of this method through case studies of the swelling behaviors of two representative systems in which swelling upon adsorption of water is of importance: single-walled aluminosilicate nanotube bundles and cesium montmorillonite. The resulting predictions show that the aluminosilicate nanotube bundles swell significantly with increasing interstitial adsorption and that the layer spacing of cesium montmorillonite expands up to about 12.5 Å, giving good agreement with experiments. The method is applicable to a wide range of flexible nanoporous materials, such as zeolites, metal-organic frameworks, and layered oxide materials, when candidate structures can be defined and a force field to describe the material is available.

Single-ensemble nonequilibrium path-sampling estimates of free energy differences.

We introduce a straightforward, single-ensemble, path sampling approach to calculate free energy differences based on Jarzynski's relation. For a two-dimensional "toy" test system, the new (minimally optimized) method performs roughly one hundred times faster than either optimized "traditional" Jarzynski calculations or conventional thermodynamic integration. The simplicity of the underlying formalism suggests the approach will find broad applicability in molecular systems.

Large Timesteps in First-Principles Molecular Dynamics Simulations

A number of first-principles molecular dynamics simulations of biologically important systems, e.g. molecular liquids and biomolecules in aqueous solution, are emerging in the past few years. Since these calculations are computationally very demanding, clever use of CPU time is crucial. One way to accelerate the simulations is to use the largest possible timestep, which is mainly limited by the highest vibrational frequency of the system. In the systems mentioned above, the highest frequencies arise from stretching motion of the covalent bonds, which are only weakly coupled to other modes in most cases. Therefore, by freezing these fast modes, we can use much larger timesteps without changing the long-time behavior of the system. Bond-angle constraints, on the other hand, will not be considered in this article, since they can significantly alter the properties of the system, with little gain in efficiency. Although the constrained molecular dynamics (CMD) is a well-established technique for classical simulations, its applications within the first-principles calculations are limited to the estimation of free-energy differences associated with chemical reactions. In these simulations, constraints were imposed only on the reaction coordinates, and the time evolution of the system was followed by the fictitious dynamics. In contrast, CMD has rarely been used for the Born-Oppenheimer dynamics (BOD), in which the electronic ground state is reached at every ionic step. In general, the number of iterations for electronic minimization increases for larger timesteps in BOD, since the quality of the extrapolated wavefunctions deteriorates. Therefore, it is not obvious if CMD provides any improvement of the performance. In this article, we investigate the efficiency of CMD within BOD as well as the effects of using deuterium atoms in place of hydrogen atoms. We have implemented CMD within our first-principles molecular dynamics code, FEMTECK, along the lines of RATTLE algorithm. Note that the computational efforts to calculate the constraint forces are negligible compared with those for the (unconstrained) ionic forces in the first-principles calculations. Therefore, we solve the nonlinear equations for the bond lengths (eqs. (A1) and (A3) of ref. 15) by the Newton-Raphson method, which, at the cost of inverting small matrices,

Free energy changes associated with amino acid substitution in proteins

The estimation of free energy differences from computer simulation of macromolecular systems is important for rational strategies for drug design and for protein engineering. As an example of one mutation, we have studied the free energy change resulting from the conversion of a polar group (OH) to an apolar group (CH3) in aqueous solution. We have estimated the effect of various local environments on the magnitude of the free energy difference and find that significant environmental effects are found. We have also studied the reliability of the results in detail.

Efficient estimation of free energy differences from Monte Carlo data

Near-optimal strategies are developed for estimating the free energy difference between two canonical ensembles, given a Metropolis-type Monte Carlo program for sampling each one. The estimation strategy depends on the extent of overlap between the two ensembles, on the smoothness of the density-of-states as a function of the difference potential, and on the relative Monte Carlo sampling costs, per statistically independent data point. The best estimate of the free energy difference is usually obtained by dividing the available computer time approximately equally between the two ensembles; its efficiency (variance x computer time)-1 is never less, and may be several orders of magnitude greater, than that obtained by sampling only one ensemble, as is done in perturbation theory.