Stratified UWHAM and Its Stochastic Approximation for Multicanonical Simulations Which Are Far from Equilibrium

Enhancing torsional sampling using fully adaptive simulated tempering

Enhanced sampling algorithms are indispensable when working with highly disconnected multimodal distributions. An important application of these is the conformational exploration of particular internal degrees of freedom of molecular systems. However, despite the existence of many commonly used enhanced sampling algorithms to explore these internal motions, they often rely on system-dependent parameters, which negatively impact efficiency and reproducibility. Here, we present fully adaptive simulated tempering (FAST), a variation of the irreversible simulated tempering algorithm, which continuously optimizes the number, parameters, and weights of intermediate distributions to achieve maximally fast traversal over a space defined by the change in a predefined thermodynamic control variable such as temperature or an alchemical smoothing parameter. This work builds on a number of previously published methods, such as sequential Monte Carlo, and introduces a novel parameter optimization procedure that can, in principle, be used in any expanded ensemble algorithms. This method is validated by being applied on a number of different molecular systems with high torsional kinetic barriers. We also consider two different soft-core potentials during the interpolation procedure and compare their performance. We conclude that FAST is a highly efficient algorithm, which improves simulation reproducibility and can be successfully used in a variety of settings with the same initial hyperparameters.

Biomolecular dynamics in the 21st century

The relevance of motions in biological macromolecules has been clear since the early structural analyses of proteins by X-ray crystallography. Computer simulations have been applied to provide a deeper understanding of the dynamics of biological macromolecules since 1976, and are now a standard tool in many labs working on the structure and function of biomolecules. In this mini-review we highlight some areas of current interest and active development for simulations, in particular all-atom molecular dynamics simulations.

Stochastic Approximation to MBAR and TRAM: Batchwise Free Energy Estimation

The dynamics of molecules are governed by rare event transitions between long-lived (metastable) states. To explore these transitions efficiently, many enhanced sampling protocols have been introduced that involve using simulations with biases or changed temperatures. Two established statistically optimal estimators for obtaining unbiased equilibrium properties from such simulations are the multistate Bennett acceptance ratio (MBAR) and the transition-based reweighting analysis method (TRAM). Both MBAR and TRAM are solved iteratively and can suffer from long convergence times. Here, we introduce stochastic approximators (SA) for both estimators, resulting in SAMBAR and SATRAM, which are shown to converge faster than their deterministic counterparts, without significant accuracy loss. Both methods are demonstrated on different molecular systems.

Hydrophobicity of arginine leads to reentrant liquid-liquid phase separation behaviors of arginine-rich proteins

Intrinsically disordered proteins rich in cationic amino acid groups can undergo Liquid-Liquid Phase Separation (LLPS) in the presence of charge-balancing anionic counterparts. Arginine and Lysine are the two most prevalent cationic amino acids in proteins that undergo LLPS, with arginine-rich proteins observed to undergo LLPS more readily than lysine-rich proteins, a feature commonly attributed to arginine's ability to form stronger cation-π interactions with aromatic groups. Here, we show that arginine's ability to promote LLPS is independent of the presence of aromatic partners, and that arginine-rich peptides, but not lysine-rich peptides, display re-entrant phase behavior at high salt concentrations. We further demonstrate that the hydrophobicity of arginine is the determining factor giving rise to the reentrant phase behavior and tunable viscoelastic properties of the dense LLPS phase. Controlling arginine-induced reentrant LLPS behavior using temperature and salt concentration opens avenues for the bioengineering of stress-triggered biological phenomena and drug delivery systems.

Iterative Landmark-Based Umbrella Sampling (ILBUS) Protocol for Sampling of Conformational Space of Biomolecules

Computer simulations of biomolecules such as molecular dynamics often suffer from insufficient sampling. Due to limited computational resources, insufficient sampling prevents obtaining proper equilibrium distributions of observed properties. To deal with this problem, we proposed a simulation protocol for efficient resampling of collected off-equilibrium trajectories. These trajectories are utilized for the initial mapping of the conformational space, which is later properly resampled by the introduced Iterative Landmark-Based Umbrella Sampling (ILBUS) method. Reconstruction of static equilibrium properties is achieved by the multistate Bennett acceptance ratio (MBAR) method, which enables efficient use of simulated data. The ILBUS protocol is geometry-based and does not demand any additional collective variable or a dimensional-reduction technique. The only requirement is a set of suitably spaced reference conformations, which serve as landmarks in the mapped conformational space. Additionally, the ILBUS protocol encompasses an iterative process that optimizes the force constant used in the umbrella sampling simulation. Such tuning is an inherent feature of the protocol and does not need to be performed by the user in advance. Furthermore, even the simulations with suboptimal force constants can be used in estimates by MBAR. We demonstrate the feasibility and the performance of this approach in the study of the conformational landscape of the alanine dipeptide, met-enkephalin, and adenylate kinase.

Protein Loop Conformational Free Energy Changes via an Alchemical Path without Reaction Coordinates

We introduce a method called restrain-free energy perturbation-release 2.0 (R-FEP-R 2.0) to estimate conformational free energy changes of protein loops via an alchemical path. R-FEP-R 2.0 is a generalization of the method called restrain-free energy perturbation-release (R-FEP-R) that can only estimate conformational free energy changes of protein side chains but not loops. The reorganization of protein loops is a central feature of many biological processes. Unlike other advanced sampling algorithms such as umbrella sampling and metadynamics, R-FEP-R and R-FEP-R 2.0 do not require predetermined collective coordinates and transition pathways that connect the two endpoint conformational states. The R-FEP-R 2.0 method was applied to estimate the conformational free energy change of a β-turn flip in the protein ubiquitin. The result obtained by R-FEP-R 2.0 agrees with the benchmarks very well. We also comment on problems commonly encountered when applying umbrella sampling to calculate protein conformational free energy changes.

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

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

The UWHAM and SWHAM Software Package

We introduce the UWHAM (binless weighted histogram analysis method) and SWHAM (stochastic UWHAM) software package that can be used to estimate the density of states and free energy differences based on the data generated by multi-state simulations. The programs used to solve the UWHAM equations are written in the C++ language and operated via the command line interface. In this paper, first we review the theoretical bases of UWHAM, its stochastic solver RE-SWHAM (replica exchange-like SWHAM)and ST-SWHAM (serial tempering-like SWHAM). Then we provide a tutorial with examples that explains how to apply the UWHAM program package to analyze the data generated by different types of multi-state simulations: umbrella sampling, replica exchange, free energy perturbation simulations, etc. The tutorial examples also show that the UWHAM equations can be solved stochastically by applying the RE-SWHAM and ST-SWHAM programs when the data ensemble is large. If the simulations at some states are far from equilibrium, the Stratified RE-SWHAM program can be applied to obtain the equilibrium distribution of the state of interest. All the source codes and the tutorial examples are available from our group's web page: https://ronlevygroup.cst.temple.edu/software/UWHAM_and_SWHAM_webpage/index.html .

The Excess Chemical Potential of Water at the Interface with a Protein from End Point Simulations

We use end point simulations to estimate the excess chemical potential of water in the homogeneous liquid and at the interface with a protein in solution. When the pure liquid is taken as the reference, the excess chemical potential of interfacial water is the difference between the solvation free energy of a water molecule at the interface and in the bulk. Using the homogeneous liquid as an example, we show that the solvation free energy for growing a water molecule can be estimated by applying UWHAM to the simulation data generated from the initial and final states (i.e., "the end points") instead of multistate free energy perturbation simulations because of the possible overlaps of the configurations sampled at the end points. Then end point simulations are used to estimate the solvation free energy of water at the interface with a protein in solution. The estimate of the solvation free energy at the interface from two simulations at the end points agrees with the benchmark using 32 states within a 95% confidence interval for most interfacial locations. The ability to accurately estimate the excess chemical potential of water from end point simulations facilitates the statistical thermodynamic analysis of diverse interfacial phenomena. Our focus is on analyzing the excess chemical potential of water at protein receptor binding sites with the goal of using this information to assist in the design of tight binding ligands.