Calculation Of Free Energy Differences Research Articles

Calculations of Absolute Solvation Free Energies with Transformato─Application to the FreeSolv Database Using the CGenFF Force Field.

We recently introduced transformato, an open-source Python package for the automated setup of large-scale calculations of relative solvation and binding free energy differences. Here, we extend the capabilities of transformato to the calculation of absolute solvation free energy differences. After careful validation against the literature results and reference calculations with the PERT module of CHARMM, we used transformato to compute absolute solvation free energies for most molecules in the FreeSolv database (621 out of 642). The force field parameters were obtained with the program cgenff (v2.5.1), which derives missing parameters from the CHARMM general force field (CGenFF v4.6). A long-range correction for the Lennard-Jones interactions was added to all computed solvation free energies. The mean absolute error compared to the experimental data is 1.12 kcal/mol. Our results allow a detailed comparison between the AMBER and CHARMM general force fields and provide a more in-depth understanding of the capabilities and limitations of the CGenFF small molecule parameters.

OneOPES, a Combined Enhanced Sampling Method to Rule Them All.

Enhanced sampling techniques have revolutionized molecular dynamics (MD) simulations, enabling the study of rare events and the calculation of free energy differences in complex systems. One of the main families of enhanced sampling techniques uses physical degrees of freedom called collective variables (CVs) to accelerate a system's dynamics and recover the original system's statistics. However, encoding all the relevant degrees of freedom in a limited number of CVs is challenging, particularly in large biophysical systems. Another category of techniques, such as parallel tempering, simulates multiple replicas of the system in parallel, without requiring CVs. However, these methods may explore less relevant high-energy portions of the phase space and become computationally expensive for large systems. To overcome the limitations of both approaches, we propose a replica exchange method called OneOPES that combines the power of multireplica simulations and CV-based enhanced sampling. This method efficiently accelerates the phase space sampling without the need for ideal CVs, extensive parameters fine tuning nor the use of a large number of replicas, as demonstrated by its successful applications to protein-ligand binding and protein folding benchmark systems. Our approach shows promise as a new direction in the development of enhanced sampling techniques for molecular dynamics simulations, providing an efficient and robust framework for the study of complex and unexplored problems.

Exploring Routes to Enhance the Calculation of Free Energy Differences via Non-Equilibrium Work SQM/MM Switching Simulations Using Hybrid Charge Intermediates between MM and SQM Levels of Theory or Non-Linear Switching Schemes.

Non-equilibrium work switching simulations and Jarzynski's equation are a reliable method for computing free energy differences, ΔAlow→high, between two levels of theory, such as a pure molecular mechanical (MM) and a quantum mechanical/molecular mechanical (QM/MM) description of a system of interest. Despite the inherent parallelism, the computational cost of this approach can quickly become very high. This is particularly true for systems where the core region, the part of the system to be described at different levels of theory, is embedded in an environment such as explicit solvent water. We find that even for relatively simple solute-water systems, switching lengths of at least 5 ps are necessary to compute ΔAlow→high reliably. In this study, we investigate two approaches towards an affordable protocol, with an emphasis on keeping the switching length well below 5 ps. Inserting a hybrid charge intermediate state with modified partial charges, which resembles the charge distribution of the desired high level, makes it possible to obtain reliable calculations with 2 ps switches. Attempts using step-wise linear switching paths, on the other hand, did not lead to improvement, i.e., a faster convergence for all systems. To understand these findings, we analyzed the solutes' properties as a function of the partial charges used and the number of water molecules in direct contact with the solute, and studied the time needed for water molecules to reorient themselves upon a change in the solute's charge distribution.

Computing Free Energies with Fluctuation Relations on Quantum Computers.

As a central thermodynamic property, free energy enables the calculation of virtually any equilibrium property of a physical system, allowing for the construction of phase diagrams and predictions about transport, chemical reactions, and biological processes. Thus, methods for efficiently computing free energies, which in general is a difficult problem, are of great interest to broad areas of physics and the natural sciences. The majority of techniques for computing free energies target classical systems, leaving the computation of free energies in quantum systems less explored. Recently developed fluctuation relations enable the computation of free energy differences in quantum systems from an ensemble of dynamic simulations. While performing such simulations is exponentially hard on classical computers, quantum computers can efficiently simulate the dynamics of quantum systems. Here, we present an algorithm utilizing a fluctuation relation known as the Jarzynski equality to approximate free energy differences of quantum systems on a quantum computer. We discuss under which conditions our approximation becomes exact, and under which conditions it serves as a strict upper bound. Furthermore, we successfully demonstrate a proof of concept of our algorithm using the transverse field Ising model on a real quantum processor. As quantum hardware continues to improve, we anticipate that our algorithm will enable computation of free energy differences for a wide range of quantum systems useful across the natural sciences.

Optimizing the Calculation of Free Energy Differences in Nonequilibrium Work SQM/MM Switching Simulations.

A key step during indirect alchemical free energy simulations using quantum mechanical/molecular mechanical (QM/MM) hybrid potential energy functions is the calculation of the free energy difference ΔAlow→high between the low level (e.g., pure MM) and the high level of theory (QM/MM). A reliable approach uses nonequilibrium work (NEW) switching simulations in combination with Jarzynski’s equation; however, it is computationally expensive. In this study, we investigate whether it is more efficient to use more shorter switches or fewer but longer switches. We compare results obtained with various protocols to reference free energy differences calculated with Crooks’ equation. The central finding is that fewer longer switches give better converged results. As few as 200 sufficiently long switches lead to ΔAlow→high values in good agreement with the reference results. This optimized protocol reduces the computational cost by a factor of 40 compared to earlier work. We also describe two tools/ways of analyzing the raw data to detect sources of poor convergence. Specifically, we find it helpful to analyze the raw data (work values from the NEW switching simulations) in a quasi-time series-like manner. Principal component analysis helps to detect cases where one or more conformational degrees of freedom are different at the low and high level of theory.

Relative free-energy calculations for scaffold hopping-type transformations with an automated RE-EDS sampling procedure

The calculation of relative free-energy differences between different compounds plays an important role in drug design to identify potent binders for a given protein target. Most rigorous methods based on molecular dynamics simulations estimate the free-energy difference between pairs of ligands. Thus, the comparison of multiple ligands requires the construction of a “state graph”, in which the compounds are connected by alchemical transformations. The computational cost can be optimized by reducing the state graph to a minimal set of transformations. However, this may require individual adaptation of the sampling strategy if a transformation process does not converge in a given simulation time. In contrast, path-free methods like replica-exchange enveloping distribution sampling (RE-EDS) allow the sampling of multiple states within a single simulation without the pre-definition of alchemical transition paths. To optimize sampling and convergence, a set of RE-EDS parameters needs to be estimated in a pre-processing step. Here, we present an automated procedure for this step that determines all required parameters, improving the robustness and ease of use of the methodology. To illustrate the performance, the relative binding free energies are calculated for a series of checkpoint kinase 1 inhibitors containing challenging transformations in ring size, opening/closing, and extension, which reflect changes observed in scaffold hopping. The simulation of such transformations with RE-EDS can be conducted with conventional force fields and, in particular, without soft bond-stretching terms.

Reweighted Jarzynski Sampling: Acceleration of Rare Events and Free Energy Calculation with a Bias Potential Learned from Nonequilibrium Work.

We introduce a simple enhanced sampling approach for the calculation of free energy differences and barriers along a one-dimensional reaction coordinate. First, a small number of short nonequilibrium simulations are carried out along the reaction coordinate, and the Jarzynski equality is used to learn an approximate free energy surface from the nonequilibrium work distribution. This free energy estimate is represented in a compact form as an artificial neural network and used as an external bias potential to accelerate rare events in a subsequent molecular dynamics simulation. The final free energy estimate is then obtained by reweighting the equilibrium probability distribution of the reaction coordinate sampled under the influence of the external bias. We apply our reweighted Jarzynski sampling recipe to four processes of varying scales and complexities─spanning chemical reaction in the gas phase, pair association in solution, and droplet nucleation in supersaturated vapor. In all cases, we find reweighted Jarzynski sampling to be a very efficient strategy, resulting in rapid convergence of the free energy to high precision.

Harnessing fluctuation theorems to discover free energy and dissipation potentials from non-equilibrium data

The Jarzynski relation, as an equality form of the second law of thermodynamics, represents an exact thermodynamic statement that is valid arbitrarily far away from equilibrium. This remarkable relation directly links the equilibrium free energy difference between two states and the probability distribution of the work done along a process that drives the system from one state to the other. Here, we leverage the Jarzynski equality and a local equilibrium assumption, to go beyond the calculation of free energy differences and also extract the dissipation potential from additional measurements of kinematic field variables (strain and velocity fields). The proposed strategy is exemplified over pulling experiments of mass–spring models obeying overdamped Langevin dynamics, which is a prototype for nanorods, fibrous macro-molecules and the Rouse model of polymers. Different interaction potentials, fluid viscosities and bath temperatures are studied, so as to intrinsically control how close or far away the system is from equilibrium. The data-inferred continuum models are then validated against processes governed by different pulling protocols, thereby demonstrating their predictive capability. The methods presented here represent a first step toward full material characterization from non-equilibrium data of macroscopic observables, which could potentially be obtained from experimental observations.

Free Energy Computation of Particles with Membrane-Mediated Interactions Via Langevin Dynamics

We apply well-established concepts of Langevin sampling to derive a new class of algorithms for the efficient computation of free energy differences of fluctuating particles embedded in a 'fast' membrane, i.e., a membrane that instantaneously adapts to varying particle positions. A geometric potential accounting for membrane-mediated particle interaction is derived in the framework of variational hybrid models for particles in membranes. Recent explicit representations of the gradient of the geometric interaction potential allows to apply well-known gradient based Markov Chain Monte-Carlo (MCDC) methods such as Langevin-based sampling.

Conformational Free-Energy Differences of Large Solvated Systems with the Focused Confinement Method.

The focused confinement method (FCM) is a reaction coordinate-free simulation approach for the calculation of conformational free-energy differences in explicit solvent. The method uses reference states for the conformations of interest, partitions the solute into conformationally active and inactive regions, and requires the calculation of desolvation free energies of mixed harmonic-anharmonic states as part of its procedure. The reference states and partitioning affect the speed of convergence of FCM's constituent simulations in opposing manners, but in the thermodynamic limit, they have no effect on calculated conformational free-energy differences. To aid fast convergence of large systems, a general procedure to quickly partition and construct reference states is introduced. With this, two sets of reference states and associated partitionings were constructed for the closed and open conformation of triosephosphate isomerase (TIM). Despite TIM's size, highly converged desolvation free energies were readily obtained from standard free-energy perturbation simulations because the mixed harmonic-anharmonic states are heavily rigidified. FCM-calculated free-energy differences for loop closing matched the experimental value for both reference sets. The insensitivity to reference states and associated partitionings favors reference states that merely reflect main structural differences for which convergence is faster. The calculations demonstrate the accuracy and robustness of FCM for large systems.

Toward Automated Free Energy Calculation with Accelerated Enveloping Distribution Sampling (A-EDS).

Free-energy perturbation (FEP) methods are commonly used in drug design to calculate relative binding free energies of different ligands to a common host protein. Alchemical ligand transformations are usually performed in multiple steps which need to be chosen carefully to ensure sufficient phase-space overlap between neighboring states. With one-step or single-step FEP techniques, a single reference state is designed that samples phase-space not only representative of a full transformation but also ideally resembles multiple ligand end states and hence allows for efficient multistate perturbations. Enveloping distribution sampling (EDS) is one example for such a method in which the reference state is created by a mathematical combination of the different ligand end states based on solid statistical mechanics. We have recently proposed a novel approach to EDS which enables efficient barrier crossing between the different end states, termed accelerated EDS (A-EDS). In this work, we further simplify the parametrization of the A-EDS reference state and demonstrate the automated calculation of multiple free-energy differences between different ligands from a single simulation in three different well-described drug design model systems.

Phase equilibrium of liquid water and hexagonal ice from enhanced sampling molecular dynamics simulations.

We study the phase equilibrium between liquid water and ice Ih modeled by the TIP4P/Ice interatomic potential using enhanced sampling molecular dynamics simulations. Our approach is based on the calculation of ice Ih-liquid free energy differences from simulations that visit reversibly both phases. The reversible interconversion is achieved by introducing a static bias potential as a function of an order parameter. The order parameter was tailored to crystallize the hexagonal diamond structure of oxygen in ice Ih. We analyze the effect of the system size on the ice Ih-liquid free energy differences, and we obtain a melting temperature of 270 K in the thermodynamic limit. This result is in agreement with estimates from thermodynamic integration (272 K) and coexistence simulations (270 K). Since the order parameter does not include information about the coordinates of the protons, the spontaneously formed solid configurations contain proton disorder as expected for ice Ih.

The Conveyor Belt Umbrella Sampling (CBUS) Scheme: Principle and Application to the Calculation of the Absolute Binding Free Energies of Alkali Cations to Crown Ethers.

We recently introduced a method called conveyor belt (CB) thermodynamic integration (TI) for the calculation of alchemical free-energy differences based on molecular dynamics simulations. In the present work, the CBTI approach is generalized to conformational free-energy changes, i.e., to the determination of the potential of mean force (PMF) along a conformational coordinate ξ of interest. The proposed conveyor belt umbrella sampling (CBUS) scheme relies on the parallel simulation of K replicas k = 0,1, ..., K - 1 of the system, with K even. For each replica k, the instantaneous value of ξ is restrained to an anchor value λk. The latter anchor points are equally spaced along a forward-turn-backward-turn path (i.e., a CB) between two extreme values defining the ξ-range of interest. The rotation of the CB is controlled by a variable Λ (range from 0 to 2π) which evolves dynamically along the simulation. The evolution of Λ results from the forces exerted by the restraining potentials on the anchor points, taken equal and opposite to those they exert on the replicas. Because these forces tend to cancel out along the CB, the dynamics of Λ is essentially diffusive, and the continuous distribution of ξ-values sampled by the replica system is automatically close to homogeneous. The latter feature represents an advantage over direct counting (DCNT) and traditional umbrella sampling (TRUS), shared to some extent with replica-exchange umbrella sampling (REUS). In this work, the CBUS scheme is introduced and compared to the three latter schemes in the calculation of 45 standard absolute binding free energies. These correspond to the binding of five alkali cations to three crown ethers in three solvents. Different free-energy estimators are considered for the PMF calculation, and the calculated values are also compared to those of a previous study relying on an alchemical path, as well as to experimental data.

Steering a solute between coexisting solvation states: Revisiting nonequilibrium work relations and the calculation of free energy differences.

By analogy with single-molecule pulling experiments, we present a computational framework to obtain free energy differences between complex solvation states. To illustrate our approach, we focus on the calculation of solvation free energies (SFEs). However, the method can be readily extended to cases involving more complex solutes and solvation conditions as well as to the calculation of binding free energies. The main idea is to drag the solute across the simulation box where atomistic and ideal gas representations of the solvent coexist at constant temperature and chemical potential. At finite pulling speeds, the resulting work allows one to extract SFEs via nonequilibrium relations, whereas at infinitely slow pulling speeds, this process becomes equivalent to the thermodynamic integration method. Results for small molecules well agree with literature data and pave the way to systematic studies of arbitrarily large and complex molecules.

Accelerating QM/MM Free Energy Computations via Intramolecular Force Matching.

The calculation of free energy differences between levels of theory has numerous potential pitfalls. Chief among them is the lack of overlap, i.e., ensembles generated at one level of theory (e.g., "low") not being good approximations of ensembles at the other (e.g., "high"). Numerous strategies have been devised to mitigate this issue. However, the most straightforward approach is to ensure that the "low" level ensemble more closely resembles that of the "high". Ideally, this is done without increasing computational cost. Herein, we demonstrate that by reparametrizing classical intramolecular potentials to reproduce high level forces (i.e., force matching) configurational overlap between a "low" (i.e., classical) and "high" (i.e., quantum) level can be significantly improved. This procedure is validated on two test cases and results in vastly improved convergence of free energy simulations.

Accelerated Enveloping Distribution Sampling: Enabling Sampling of Multiple End States while Preserving Local Energy Minima.

Enveloping distribution sampling (EDS) is an efficient approach to calculate multiple free-energy differences from a single molecular dynamics (MD) simulation. However, the construction of an appropriate reference-state Hamiltonian that samples all states efficiently is not straightforward. We propose a novel approach for the construction of the EDS reference-state Hamiltonian, related to a previously described procedure to smoothen energy landscapes. In contrast to previously suggested EDS approaches, our reference-state Hamiltonian preserves local energy minima of the combined end-states. Moreover, we propose an intuitive, robust and efficient parameter optimization scheme to tune EDS Hamiltonian parameters. We demonstrate the proposed method with established and novel test systems and conclude that our approach allows for the automated calculation of multiple free-energy differences from a single simulation. Accelerated EDS promises to be a robust and user-friendly method to compute free-energy differences based on solid statistical mechanics.

Unsupervised Calculation of Free Energy Barriers in Large Crystalline Systems.

The calculation of free energy differences for thermally activated mechanisms in the solid state are routinely hindered by the inability to define a set of collective variable functions that accurately describe the mechanism under study. Even when possible, the requirement of descriptors for each mechanism under study prevents implementation of free energy calculations in the growing range of automated material simulation schemes. We provide a solution, deriving a path-based, exact expression for free energy differences in the solid state which does not require a converged reaction pathway, collective variable functions, Gram matrix evaluations, or probability flux-based estimators. The generality and efficiency of our method is demonstrated on a complex transformation of C15 interstitial defects in iron and double kink nucleation on a screw dislocation in tungsten, the latter system consisting of more than 120 000 atoms. Both cases exhibit significant anharmonicity under experimentally relevant temperatures.

Calculation of conformational free energies by confinement simulations in explicit water with implicit desolvation

The confinement method is a robust and conceptually simple free energy simulation method that allows the calculation of conformational free energy differences between highly dissimilar states. Application of the method to explicitly solvated systems requires a multi-stage simulation protocol for the calculation of desolvation free energies. Here we show that these desolvation free energies can be readily obtained from an implicit treatment, which is simpler and less costly. The accuracy and robustness of this protocol was shown by the calculation of conformational free energy differences of a series of explicitly solvated test systems. Given the accuracy and ease by which these free energy differences were obtained, the confinement method is promising for the treatment of conformational changes in large and complex systems.

Using force-based adaptive resolution simulations to calculate solvation free energies of amino acid sidechain analogues.

The calculation of free energy differences is a crucial step in the characterization and understanding of the physical properties of biological molecules. In the development of efficient methods to compute these quantities, a promising strategy is that of employing a dual-resolution representation of the solvent, specifically using an accurate model in the proximity of a molecule of interest and a simplified description elsewhere. One such concurrent multi-resolution simulation method is the Adaptive Resolution Scheme (AdResS), in which particles smoothly change their resolution on-the-fly as they move between different subregions. Before using this approach in the context of free energy calculations, however, it is necessary to make sure that the dual-resolution treatment of the solvent does not cause undesired effects on the computed quantities. Here, we show how AdResS can be used to calculate solvation free energies of small polar solutes using Thermodynamic Integration (TI). We discuss how the potential-energy-based TI approach combines with the force-based AdResS methodology, in which no global Hamiltonian is defined. The AdResS free energy values agree with those calculated from fully atomistic simulations to within a fraction of kBT. This is true even for small atomistic regions whose size is on the order of the correlation length, or when the properties of the coarse-grained region are extremely different from those of the atomistic region. These accurate free energy calculations are possible because AdResS allows the sampling of solvation shell configurations which are equivalent to those of fully atomistic simulations. The results of the present work thus demonstrate the viability of the use of adaptive resolution simulation methods to perform free energy calculations and pave the way for large-scale applications where a substantial computational gain can be attained.

Balancing Accuracy and Cost of Confinement Simulations by Interpolation and Extrapolation of Confinement Energies.

Improvements to the confinement method for the calculation of conformational free energy differences are presented. By taking advantage of phase space overlap between simulations at different frequencies, significant gains in accuracy and speed are reached. The optimal frequency spacing for the simulations is obtained from extrapolations of the confinement energy, and relaxation time analysis is used to determine time steps, simulation lengths, and friction coefficients. At postprocessing, interpolation of confinement energies is used to significantly reduce discretization errors in the calculation of conformational free energies. The efficiency of this protocol is illustrated by applications to alanine n-peptides and lactoferricin. For the alanine-n-peptide, errors were reduced between 2- and 10-fold and sampling times between 8- and 67-fold, while for lactoferricin the long sampling times at low frequencies were reduced 10-100-fold.