A hybrid recursion method to robustly ensure convergence efficiencies in the simulated scaling based free energy simulations

An introduction to best practices in free energy calculations

Free energy calculations are extremely useful for investigating small-molecule biophysical properties such as protein-ligand binding affinities and partition coefficients. However, these calculations are also notoriously difficult to implement correctly. In this chapter, we review standard methods for computing free energy via simulation, discussing current best practices and examining potential pitfalls for computational researchers performing them for the first time. We include a variety of examples and tips for how to set up and conduct these calculations, including applications to relative binding affinities and small-molecule solvation free energies.

Generalized essential energy space random walks to more effectively accelerate solute sampling in aqueous environment

Molecular dynamics sampling can be enhanced via the promoting of potential energy fluctuations, for instance, based on a Hamiltonian modified with the addition of a potential-energy-dependent biasing term. To overcome the diffusion sampling issue, which reveals the fact that enlargement of event-irrelevant energy fluctuations may abolish sampling efficiency, the essential energy space random walk (EESRW) approach was proposed earlier. To more effectively accelerate the sampling of solute conformations in aqueous environment, in the current work, we generalized the EESRW method to a two-dimension-EESRW (2D-EESRW) strategy. Specifically, the essential internal energy component of a focused region and the essential interaction energy component between the focused region and the environmental region are employed to define the two-dimensional essential energy space. This proposal is motivated by the general observation that in different conformational events, the two essential energy components have distinctive interplays. Model studies on the alanine dipeptide and the aspartate-arginine peptide demonstrate sampling improvement over the original one-dimension-EESRW strategy; with the same biasing level, the present generalization allows more effective acceleration of the sampling of conformational transitions in aqueous solution. The 2D-EESRW generalization is readily extended to higher dimension schemes and employed in more advanced enhanced-sampling schemes, such as the recent orthogonal space random walk method.

Best practices in free energy calculations for drug design

Free energy calculations are increasingly of interest for computing biophysical properties of novel small molecules of interest in drug design, such as protein-ligand binding affinities and small molecule partition coefficients. However, these calculations are also notoriously difficult to implement correctly. In this article, we review standard methods for computing free energy differences via simulation, discuss current best practices, and examine potential pitfalls for computational researchers without extensive experience in such calculations. We include a variety of examples and tips for how to set up and conduct these calculations, including applications to relative binding affinities and absolute binding free energies.

Force-momentum-based self-guided Langevin dynamics: a rapid sampling method that approaches the canonical ensemble

The self-guided Langevin dynamics (SGLD) is a method to accelerate conformational searching. This method is unique in the way that it selectively enhances and suppresses molecular motions based on their frequency to accelerate conformational searching without modifying energy surfaces or raising temperatures. It has been applied to studies of many long time scale events, such as protein folding. Recent progress in the understanding of the conformational distribution in SGLD simulations makes SGLD also an accurate method for quantitative studies. The SGLD partition function provides a way to convert the SGLD conformational distribution to the canonical ensemble distribution and to calculate ensemble average properties through reweighting. Based on the SGLD partition function, this work presents a force-momentum-based self-guided Langevin dynamics (SGLDfp) simulation method to directly sample the canonical ensemble. This method includes interaction forces in its guiding force to compensate the perturbation caused by the momentum-based guiding force so that it can approximately sample the canonical ensemble. Using several example systems, we demonstrate that SGLDfp simulations can approximately maintain the canonical ensemble distribution and significantly accelerate conformational searching. With optimal parameters, SGLDfp and SGLD simulations can cross energy barriers of more than 15 kT and 20 kT, respectively, at similar rates for LD simulations to cross energy barriers of 10 kT. The SGLDfp method is size extensive and works well for large systems. For studies where preserving accessible conformational space is critical, such as free energy calculations and protein folding studies, SGLDfp is an efficient approach to search and sample the conformational space.

Developing hybrid approaches to predict pKa values of ionizable groups

Accurate predictions of pKa values of titratable groups require taking into account all relevant processes associated with the ionization/deionization. Frequently, however, the ionization does not involve significant structural changes and the dominating effects are purely electrostatic in origin allowing accurate predictions to be made based on the electrostatic energy difference between ionized and neutral forms alone using a static structure and the subtle structural changes be accounted by using dielectric constant larger than two. On another hand, if the change of the charge state is accompanied by a large structural reorganization of the target protein, then the relevant conformational changes have to be explicitly taken into account in the pKa calculations. Here we report a hybrid approach that first predicts the titratable groups whose ionization is expected to cause large conformational changes, termed "problematic" residues, and then applies a special protocol on them, while the rest of the pK(a)s are predicted with rigid backbone approach as implemented in multi-conformation continuum electrostatics (MCCE) method. The backbone representative conformations for "problematic" groups are generated with either molecular dynamics simulations with charged and uncharged amino acid or with ab-initio local segment modeling. The corresponding ensembles are then used to calculate the titration curves of the "problematic" residues and then the results are averaged to obtain the corresponding pKa.

Computing Relative Free Energies of Solvation using Single Reference Thermodynamic Integration Augmented with Hamiltonian Replica Exchange

This paper introduces an efficient single-topology variant of Thermodynamic Integration (TI) for computing relative transformation free energies in a series of molecules with respect to a single reference state. The presented TI variant that we refer to as Single-Reference TI (SR-TI) combines well-established molecular simulation methodologies into a practical computational tool. Augmented with Hamiltonian Replica Exchange (HREX), the SR-TI variant can deliver enhanced sampling in select degrees of freedom. The utility of the SR-TI variant is demonstrated in calculations of relative solvation free energies for a series of benzene derivatives with increasing complexity. Noteworthy, the SR-TI variant with the HREX option provides converged results in a challenging case of an amide molecule with a high (13-15 kcal/mol) barrier for internal cis/trans interconversion using simulation times of only 1 to 4 ns.

Simultaneous escaping of explicit and hidden free energy barriers: application of the orthogonal space random walk strategy in generalized ensemble based conformational sampling

To overcome the pseudoergodicity problem, conformational sampling can be accelerated via generalized ensemble methods, e.g., through the realization of random walks along prechosen collective variables, such as spatial order parameters, energy scaling parameters, or even system temperatures or pressures, etc. As usually observed, in generalized ensemble simulations, hidden barriers are likely to exist in the space perpendicular to the collective variable direction and these residual free energy barriers could greatly abolish the sampling efficiency. This sampling issue is particularly severe when the collective variable is defined in a low-dimension subset of the target system; then the "Hamiltonian lagging" problem, which reveals the fact that necessary structural relaxation falls behind the move of the collective variable, may be likely to occur. To overcome this problem in equilibrium conformational sampling, we adopted the orthogonal space random walk (OSRW) strategy, which was originally developed in the context of free energy simulation [L. Zheng, M. Chen, and W. Yang, Proc. Natl. Acad. Sci. U.S.A. 105, 20227 (2008)]. Thereby, generalized ensemble simulations can simultaneously escape both the explicit barriers along the collective variable direction and the hidden barriers that are strongly coupled with the collective variable move. As demonstrated in our model studies, the present OSRW based generalized ensemble treatments show improved sampling capability over the corresponding classical generalized ensemble treatments.

Enhanced sampling in generalized ensemble with large gap of sampling parameter: case study in temperature space random walk

We present an efficient sampling method for computing a partition function and accelerating configuration sampling. The method performs a random walk in the lambda space, with lambda being any thermodynamic variable that characterizes a canonical ensemble such as the reciprocal temperature beta or any variable that the Hamiltonian depends on. The partition function is determined by minimizing the difference of the thermal conjugates of lambda (the energy in the case of lambda = beta), defined as the difference between the value from the dynamically updated derivatives of the partition function and the value directly measured from simulation. Higher-order derivatives of the partition function are included to enhance the Brownian motion in the lambda space. The method is much less sensitive to the system size, and to the size of lambda window than other methods. On the two dimensional Ising model, it is shown that the method asymptotically converges the partition function, and the error of the logarithm of the partition function is much smaller than the algorithm using the Wang-Landau recursive scheme. The method is also applied to off-lattice model proteins, the AB models, in which cases many low energy states are found in different models.

Random walk in orthogonal space to achieve efficient free-energy simulation of complex systems

In the past few decades, many ingenious efforts have been made in the development of free-energy simulation methods. Because complex systems often undergo nontrivial structural transition during state switching, achieving efficient free-energy calculation can be challenging. As identified earlier, the "Hamiltonian" lagging, which reveals the fact that necessary structural relaxation falls behind the order parameter move, has been a primary problem for generally low free-energy simulation efficiency. Here, we propose an algorithm by achieving a random walk in both the order parameter space and its generalized force space; thereby, the order parameter move and the required conformational relaxation can be efficiently synchronized. As demonstrated in both the alchemical transition and the conformational transition, a leapfrog improvement in free-energy simulation efficiency can be obtained; for instance, (i) it allows us to solve a notoriously challenging problem: accurately predicting the pK(a) value of a buried titratable residue, Asp-66, in the interior of the V66E staphylococcal nuclease mutant, and (ii) it allows us to gain superior efficiency over the metadynamics algorithm.