Recent Advances in Maximum Entropy Biasing Techniques for Molecular   Dynamics

Dilnoza B. Amirkulova; Andrew D. White

arXiv:1902.02252·physics.chem-ph·May 15, 2019

Recent Advances in Maximum Entropy Biasing Techniques for Molecular Dynamics

Dilnoza B. Amirkulova, Andrew D. White

PDF

TL;DR

This review discusses recent maximum entropy biasing techniques in molecular dynamics that incorporate experimental data to enhance simulation accuracy, especially for approximate force fields and multiscale processes.

Contribution

It introduces and reviews the latest experiment directed simulation (EDS) and experiment directed metadynamics (EDM) methods for improved data integration in molecular simulations.

Findings

01

Applications to ab initio water simulations

02

Incorporation of environmental fluctuations in protein complexes

03

Enhancement of RNA force fields

Abstract

This review describes recent advances by the authors and others on the topic of incorporating experimental data into molecular simulations through maximum entropy methods. Methods which incorporate experimental data improve accuracy in molecular simulation by minimally modifying the thermodynamic ensemble. This is especially important where force fields are approximate, such as when employing coarse-grain models, or where high accuracy is required, such as when attempting to mimic a multiscale self-assembly process. The authors review here the experiment directed simulation (EDS) and experiment directed metadynamics (EDM) methods that allow matching averages and distributions in simulations, respectively. Important system-specific considerations are discussed such as using enhanced sampling simultaneously, the role of pressure, treating uncertainty, and implementations of these methods.…

Equations46

Δ S_{rel} = \int d r P (r) ln \frac{P ( r )}{P ^{'} ( r )}

Δ S_{rel} = \int d r P (r) ln \frac{P ( r )}{P ^{'} ( r )}

⟨ s (r) ⟩ = \int d r P^{'} (r) s (r) = \overset{s}{^}

⟨ s (r) ⟩ = \int d r P^{'} (r) s (r) = \overset{s}{^}

L [λ, P (r), P^{'} (r)] = Δ S_{rel} - λ (\int d r P^{'} (r) s (r) - \overset{s}{^})

L [λ, P (r), P^{'} (r)] = Δ S_{rel} - λ (\int d r P^{'} (r) s (r) - \overset{s}{^})

P^{'} (r) = \frac{1}{Z ^{'}} e^{- β (U (r) + λ s (r))}

P^{'} (r) = \frac{1}{Z ^{'}} e^{- β (U (r) + λ s (r))}

U^{'} (r) = U (r) + i \sum λ_{i} s_{i} (r) + j \sum μ_{j} [v (r)]

U^{'} (r) = U (r) + i \sum λ_{i} s_{i} (r) + j \sum μ_{j} [v (r)]

\int d r δ (v (r) - v^{'}) P^{'} (r) = q (v^{'})

\int d r δ (v (r) - v^{'}) P^{'} (r) = q (v^{'})

α_{τ + 1} = α_{τ} + η_{τ} g_{τ}

α_{τ + 1} = α_{τ} + η_{τ} g_{τ}

g_{τ} = \frac{2 β}{w} (⟨ s ⟩_{τ} - \overset{s}{^}) (⟨ s^{2} ⟩_{τ} - ⟨ s ⟩_{τ}^{2})

g_{τ} = \frac{2 β}{w} (⟨ s ⟩_{τ} - \overset{s}{^}) (⟨ s^{2} ⟩_{τ} - ⟨ s ⟩_{τ}^{2})

η_{τ} = \frac{A}{\sum _{i}^{τ} g _{i}^{2}}

η_{τ} = \frac{A}{\sum _{i}^{τ} g _{i}^{2}}

μ (v)_{τ + 1} = μ (v)_{τ} + \frac{1}{q ( v _{τ} )} exp (- θ (v_{τ}, v)_{τ}) G (v, v_{τ})

μ (v)_{τ + 1} = μ (v)_{τ} + \frac{1}{q ( v _{τ} )} exp (- θ (v_{τ}, v)_{τ}) G (v, v_{τ})

\theta(v_{\tau},v)_{\tau}=\left\{\begin{array}[]{lr}1,&\textrm{non-convergent}\\ \hat{\mu}_{\tau},&q(v)\\ \beta\Delta T\mu(v_{\tau})_{\tau}/T,&q(v)^{\Delta T/(\Delta T+T)}P(v)^{T/(\Delta T+T)}\end{array}\right.

\theta(v_{\tau},v)_{\tau}=\left\{\begin{array}[]{lr}1,&\textrm{non-convergent}\\ \hat{\mu}_{\tau},&q(v)\\ \beta\Delta T\mu(v_{\tau})_{\tau}/T,&q(v)^{\Delta T/(\Delta T+T)}P(v)^{T/(\Delta T+T)}\end{array}\right.

\int d r P^{''} (r, ϵ) [s (r) + ϵ] = \overset{s}{^}

\int d r P^{''} (r, ϵ) [s (r) + ϵ] = \overset{s}{^}

g_{τ} = \frac{2 β}{w} (⟨ s ⟩_{τ} + ξ (α_{τ}) - \overset{s}{^}) (⟨ s^{2} ⟩_{τ} + \frac{\partial ξ ( α _{τ} )}{\partial α} - ⟨ s ⟩_{τ}^{2})

g_{τ} = \frac{2 β}{w} (⟨ s ⟩_{τ} + ξ (α_{τ}) - \overset{s}{^}) (⟨ s^{2} ⟩_{τ} + \frac{\partial ξ ( α _{τ} )}{\partial α} - ⟨ s ⟩_{τ}^{2})

ξ (α_{τ}) = \frac{\int d ϵ P _{0} ( ϵ ) e ^{- α_{τ} ϵ} ϵ}{\int d ϵ P _{0} ( ϵ ) e ^{- α_{τ} ϵ}}

ξ (α_{τ}) = \frac{\int d ϵ P _{0} ( ϵ ) e ^{- α_{τ} ϵ} ϵ}{\int d ϵ P _{0} ( ϵ ) e ^{- α_{τ} ϵ}}

g_{τ} = \frac{2 β}{w} (⟨ s ⟩_{τ} - \overset{s}{^}) (⟨ s^{2} ⟩_{τ} - ⟨ s ⟩_{τ}^{2}) + 2 ν Δ p_{τ} \frac{\partial Δ p _{τ}}{\partial α _{τ}}

g_{τ} = \frac{2 β}{w} (⟨ s ⟩_{τ} - \overset{s}{^}) (⟨ s^{2} ⟩_{τ} - ⟨ s ⟩_{τ}^{2}) + 2 ν Δ p_{τ} \frac{\partial Δ p _{τ}}{\partial α _{τ}}

Δ p_{τ} = - \frac{α _{τ}}{w} \frac{d s}{d V} = - \frac{α _{τ}}{w} ij \sum (\frac{\partial s}{\partial v _{ij}} \frac{\partial v _{ij}}{\partial V})

Δ p_{τ} = - \frac{α _{τ}}{w} \frac{d s}{d V} = - \frac{α _{τ}}{w} ij \sum (\frac{\partial s}{\partial v _{ij}} \frac{\partial v _{ij}}{\partial V})

ν = ν^{*} \frac{w ^{2} β ^{2}}{ρ ^{2}}

ν = ν^{*} \frac{w ^{2} β ^{2}}{ρ ^{2}}

\int d r P^{'} (r) = 1

\int d r P^{'} (r) = 1

V [P^{'} (r)] = \frac{1}{3 V} \int d r r \cdot \frac{\partial ln P ^{'} ( r )}{\partial r} P^{'} (r) = \int d r r \cdot \nabla P^{'} (r) = p_{e}

V [P^{'} (r)] = \frac{1}{3 V} \int d r r \cdot \frac{\partial ln P ^{'} ( r )}{\partial r} P^{'} (r) = \int d r r \cdot \nabla P^{'} (r) = p_{e}

\int P (r) ln \frac{P ^{'} ( r )}{P ( r )}

\int P (r) ln \frac{P ^{'} ( r )}{P ( r )}

\frac{Δ p [ P ^{'} ( r )]}{δ P ^{'} ( r )} = - \nabla [r_{1}, r_{2}, ...] = - N

\frac{Δ p [ P ^{'} ( r )]}{δ P ^{'} ( r )} = - \nabla [r_{1}, r_{2}, ...] = - N

\int d r r \cdot \nabla [P (r) + ϵ ϕ (r)] - \int d r r \cdot \nabla P (r) \approx - N \int d r ϵ ϕ (r)

\int d r r \cdot \nabla [P (r) + ϵ ϕ (r)] - \int d r r \cdot \nabla P (r) \approx - N \int d r ϵ ϕ (r)

\frac{δ L}{δ P ^{'} ( r )} = ln P^{'} (r) - ln P (r) + λ_{1} - λ_{2} N

\frac{δ L}{δ P ^{'} ( r )} = ln P^{'} (r) - ln P (r) + λ_{1} - λ_{2} N

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Recent Advances in Maximum Entropy Biasing Techniques for Molecular Dynamics

\nameD. Amirkulovaa and A. D Whitea CONTACT A. D. White. Email: [email protected] a University of Rochester, Department of Chemical Engineering, Rochester, NY, USA

Abstract

This review describes recent advances by the authors and others on the topic of incorporating experimental data into molecular simulations through maximum entropy methods. Methods which incorporate experimental data improve accuracy in molecular simulation by minimally modifying the thermodynamic ensemble. This is especially important where force fields are approximate, such as when employing coarse-grain models, or where high accuracy is required, such as when attempting to mimic a multiscale self-assembly process. The authors review here the experiment directed simulation (EDS) and experiment directed metadynamics (EDM) methods that allow matching averages and distributions in simulations, respectively. Important system-specific considerations are discussed such as using enhanced sampling simultaneously, the role of pressure, treating uncertainty, and implementations of these methods. Recent examples of EDS and EDM are reviewed including applications to ab initio molecular dynamics of water, incorporating environmental fluctuations inside of a macromolecular protein complex, improving RNA force fields, and the combination of enhanced sampling with minimal biasing to model peptides.

keywords:

molecular dynamics; maximum entropy

††articletype: Review

1 Introduction

A common task in molecular simulation is ensuring that observables of the simulation match experimentally measured values. For example, in simulations of protein structure[1] or liquids[2], quantitative agreement with experiments is the standard for assessing correctness of a model. When there is no quantitative agreement, changing the potential energy function or adding additional components to the simulation are possible ways to improve the fit. Making such changes can be an ambiguous and challenging process, especially if the potential energy function has multiple terms that can be modified. Minimal biasing techniques are a class of methods that modify a potential energy function to improve quantitative agreement with experimental values while minimizing the the change in the potential energy function. The definition of “minimal” and the way the potential energy function is modified vary from method to method.

Recent reviews of minimal biasing methods can be found in Sormanni et al. [3], Bonomi et al. [4], and Olsson et al. [5]. Table 1 in Bonomi et al. [4] provides an overview of 28 minimal biasing methods, categorizing them by whether they maximize entropy, maximize parsimony, or use Bayesian inference. These three categories correspond broadly to the criteria used to ensure that the biasing function introduces a minimal change to the potential energy function. This review focuses on two techniques developed by the authors that are categorized as maximum entropy methods: experiment directed simulation[6] (EDS) and experiment directed metadynamics[7] (EDM). EDS is for matching ensemble average scalars and EDM is for matching free energy surfaces (probability distributions of observables).

EDS, like other minimal biasing methods, modifies a potential energy function to change ensemble averages of observables to match a specific value. These observables are typically equivalent to collective variables, but are intended to be only those that area experimentally verifiable quantities. What separates EDS from other methods is that it does not use replicas and can be used to construct a continuous NVE trajectory. For example, the Bayesian landscape tilting method of Beauchamp et al. [8] relies on post-processing so that there cannot be a continuous NVE trajectory. Another example is the replica method of Lindorff-Larsen et al. [9] which relies on replica-exchange of biasing forces and thus cannot result in a continuous trajectory. This does not mean EDS is “better”; indeed these two methods appear to provide better sampling and scaling than EDS. However, the ability to compute an NVE trajectory allows dynamic observables like hydrogen-bonding lifetimes to be computed. The key results from EDS have been to improve thermodynamic observables and indirectly improve dynamic observables. For example, EDS was recently used to create a state-of-the-art DFT water model that gives near perfect agreement with X-ray scattering results, water diffusivity, and proton-hopping behavior by improving only the water oxygen-oxygen coordination number[10].

EDM is a maximum entropy method that matches ensemble probability distributions, or equivalently free energy surfaces, using prescribed functions. So far, EDM has proved most useful for matching radial distribution functions (RDFs), e.g. to then do coarse-grained modeling. EDM is less often used because it is rare to have experimental data that gives a probability distribution, other than a normal distribution (which is better treated with methods like in Hummer and Köfinger [11]). Other groups have since arrived at the same approach as EDM and typically it is now called “targeted metadynamics” because it uses the method of metadynamics to arrive at a target free energy surface[12, 13]. EDM is also equivalent to variationally enhanced sampling metadynamics[14], although there the target is a tool to improve sampling, and the biased simulation is not the goal. The name EDM was chosen by the authors because at the time there were new metadynamics methods that could target a collective variable domain[15]. Now the term “targeted” seems to apply exclusively to target distributions, and EDM is probably best described as falling under the umbrella of targeted metadynamics methods.

Here we review the maximum entropy derivations that underpin both EDS and EDM, describe some of the recent research benefiting from these methods, and discuss implementation details such as understanding the effect of EDS and EDM on the virial and accounting for uncertainty in the experimental data.

2 Theory

EDS minimally modifies an ensemble so that the ensemble average of some scalar matches a desired value, such as a value obtained from an experiment. The EDS bias is minimal because the resulting biased ensemble maximizes entropy[16, 17]. Maximum entropy and minimum relative entropy are equivalent approaches to derive these minimal biasing equations[18]. We will use the minimum relative entropy approach because it has an intuitive interpretation as a distance metric. Figure 1 illustrates how the biased ensemble is one of many choices that matches our constraints, but is as close as possible to the unbiased ensemble as measured via relative entropy.

Consider an unbiased potential energy function, $U(\vec{r})$ , which has a probability distribution of $P(\vec{r})$ under the NVT ensemble, following the Boltzmann distribution: $P(\vec{r})\propto\exp\left(-\beta U(\vec{r})\right)$ . $\vec{r}$ is the set of coordinate vectors of the $N$ particles of the ensemble and $\beta=\frac{1}{kT}$ . $k$ is the Boltzmann constant and $T$ is absolute temperature. We would like to find a biased ensemble $P^{\prime}(\vec{r})$ which is as similar as possible to the unbiased ensemble. Similarity can be defined via the relative entropy:

[TABLE]

where $P^{\prime}(\vec{r})$ is the biased ensemble probability distribution, and the integral is taken over all coordinates. Having a lower $\Delta S_{\textrm{rel}}$ means that the biased and unbiased ensembles are more similar. The biased ensemble should also have an observable average that matches a target value, which could be obtained from an experiment. This constraint is represented as

[TABLE]

where $s(\vec{r})$ is an instantaneous value for our observable (collective variable) which we are matching to $\hat{s}$ , the desired scalar value. This is sometimes called a forward model and depends only on positions. We will relax this assumption below.

These equations represent (1) a constraint (to have our average simulation value match the target $\hat{s}$ ), and (2) a scalar to minimize (the relative entropy). When these conditions are present, we can find the optimal value via the method of Lagrange multipliers. The Lagrangian is:

[TABLE]

Equation 3 can be minimized by taking the functional derivative equal to zero: $\frac{\delta\mathcal{L}}{\delta P^{\prime}(\vec{r})}$ = 0. Solving for $P^{\prime}(\vec{r})$ gives this expression for the biased ensemble:

[TABLE]

where $Z^{\prime}$ is a normalization constant. This gives an expression for the biased potential energy as $U^{\prime}(\vec{r})=U(\vec{r})+\lambda s(\vec{r})$ . This result shows that if the bias is linear in the instantaneous observable, the ensemble is minimally biased. This derivation can be repeated for multiple dimensions[16] and for functions instead of scalars[7]. The general result is:

[TABLE]

where $i$ is the index of ensemble averages that are matched to a desired value and $j$ is the index of free energy surfaces that are matched to desired functions. $\mu_{j}\left[v(\vec{r})\right]$ is a bias added that depends on $v(\vec{r})$ , another instantaneous observable, and causes the biased ensemble to match a specific distribution in the free energy surfaces $F_{t}\left[v(\vec{r})\right]$ . Specifically,

[TABLE]

where $\delta$ is the Dirac delta function and $q(v^{\prime})$ is a desired probability distribution for $v(\vec{r})$ (i.e., $q[v(\vec{r})]=-\frac{1}{\beta}\ln F_{t}\left[v(\vec{r})\right]$ ). $q(v)$ could be obtained for example from a scattering or FRET experiment[19].

The derivation above gives the form of the biased potential energy that is minimally biased, but it does not enable calculation of the Lagrange multipliers. That is the purpose of the EDS and EDM methods. EDS and EDM are time-dependent methods that change the potential energy of a simulation to arrive at the bias in Equation 5. EDS is intended to be used in a two-step process: finding the Lagrange multiplier (adaptive) and then running a standard MD simulation using the modified force field given by Equation 5.

During the adaptive phase of EDS, the Lagrange multipliers are called coupling constants to distinguish from the time-independent Lagrange multipliers. These are indicated as $\alpha_{\tau}$ where $\tau$ is a discrete step index. The Lagrange multipliers are set to be the average of $\alpha_{\tau}$ . $\alpha_{\tau}$ is defined as

[TABLE]

where $w$ is an arbitrary constant used to ensure unit homogeneity, $\left<\cdot\right>_{\tau}$ is the ensemble average between step $\tau-1$ and $\tau$ , and $\eta_{t}$ is

[TABLE]

where $A$ is a user-defined constant that controls the size of the first step. The point of $\eta_{\tau}$ is to reduce the size of the steps over time. Note that because $g_{\tau}$ is in the sum, $\left|\alpha_{0}\right|=A$ . Typically the gradients should be clipped for stability and $A$ is a natural upper/lower bound for clipping it. Thus a user of this method must choose $A$ and the time between updates. This update procedure is derived in White and Voth [6] and is based on per-coordinate infinite horizon stochastic gradient descent[20]. Hocky et al. [21] recently found that Equation 8 can be modified to use covariance in multiple dimensions to improve convergenece and that replacing Equation 7 with Levenberg–Marquardt optimization further improves convergence.

The EDM method finds the $\mu\left[v(\vec{r})\right]$ bias function in Equation 5. For compactness, we will now omit the dependence on $\vec{r}$ . Unlike EDS, EDM consists of a single phase where $\mu(v)$ changes less over time. The update equation is

[TABLE]

where $v_{\tau}$ is the value of $v(\vec{r})$ at time $\tau$ , $G$ is a kernel function (e.g., Gaussian), $q(v_{\tau})$ is the probability at position $v_{\tau}$ from the target probability distribution and $\theta(v_{\tau},v)_{\tau}$ is a function which controls convergence. Following the above update rule causes the simulation to converge to the following distributions[7, 22], depending on the choice of $\theta$ :

[TABLE]

where $\Delta T$ is the “tempering factor”[7], $P(v)$ is the marginal unbiased distribution for $v(\vec{r})$ , and $\hat{\mu}$ is the total or average of all previous $\mu_{\tau}$ . The first condition, $\theta=1$ , is similar to normal metadynamics; it reaches a distribution similar to $q(v)$ but may oscillate due to a lack of dampening in the update size. Condition two, called globally tempered, converges correctly to $q(v)$ . The last condition, locally tempered, converges to an adjustable mixture of the unbiased and target distribution. EDM is traditionally implemented as globally tempered, so that the the final distribution is indeed the target. The first condition is good for tuning parameters, since it makes progress more quickly. The final condition, locally tempered, allows a pseudo-Bayesian tuning of prior belief in the unbiased ensemble. It is pseudo-Bayesian because the ratio of influence from the prior belief to evidence is computed, not set, in Bayesian modeling.

An important consideration of both EDS and EDM is that they add potential energy during the update step which quickly becomes kinetic energy. Thus, it is important that the thermostat used to maintain constant temperature in the NVT can dissipate energy faster than it is added by the update steps. This is necessary during the adaptive phase of EDS and EDM prior to convergence.

2.1 Treating Uncertainty in Experimental Data

One complication of minimal biasing methods is uncertainty in the experimental data. For example, the experimental data could be the radius of gyration of a polymer with a reported uncertainty in the mean of 5 nm. Should the radius of gyration be matched exactly, or only to within 5 nm? It is possible to stop the adaptive phase of EDS early, for example when the average is within the uncertainty of the experimental data. An early stop is ad-hoc and part of the maximum entropy derivation unlike other methods which are built to address uncertainty[11, 23].

Cesari et al. [24] proposed a modification to the maximum entropy derivation above to address experimental uncertainty. Instead of the constraint in Equation 2, this constraint is used

[TABLE]

where $\epsilon$ is an auxiliary variable that allows deviations in the average of $s(\vec{r})$ and $P^{\prime\prime}(\cdot)$ is the probability distribution that will be solved after maximizing entropy. $\epsilon$ is a random variable from a prior distribution $P_{0}(\epsilon)$ that describes the uncertainty in the experimental data. For example, $P_{0}(\epsilon)$ could be a normal distribution or a Laplace distribution. Cesari et al. [24] show that $P^{\prime\prime}(\vec{r},\epsilon)=1/Z^{\prime\prime}P^{\prime}(\vec{r})P_{0}(\epsilon)e^{\lambda\epsilon}$ . This leads to a different update step of

[TABLE]

where $\xi(\alpha_{\tau})$ is the analytic posterior of $\epsilon$ :

[TABLE]

Returning to the radius of gyration example above, if the uncertainty ( $P_{0}(\epsilon)$ is assumed to be a normal distribution normal with $\sigma=5$ nm, then $\xi(\alpha_{\tau})=-25\alpha_{\tau}$ . This new update-step can be used to rigorously include experimental uncertainty in to the EDS method.

EDM is able to tune the relative importance of the target distribution and the unbiased ensemble with the locally-tempered variant in Equation 11. This could be used to match intuition. For example, if you believe that the target distribution which comes from experimental data is twice as accurate as the molecular dynamics simulation, you could choose $\Delta T/(\Delta T+T)=0.66$ giving about twice as much weight to the target distribution.

3 Applications of EDS and EDM

A model 1-D system is shown in Figure 2 as a probability distribution function, $P(r)$ . EDS is being used to modify the average value of $r$ to match a new set point, $\hat{r}$ . EDS adds a linear bias, whose strength is indicated with the red dot, to create the biased PDF $P^{\prime}(r)$ according to Equation 5. Notice how the features of $P(r)$ are mostly maintained in $P^{\prime}(r)$ . The bottom plot shows how each value of $\hat{r}$ corresponds to a unique biasing strength.

A more sophisticated system which demonstrates the capabilities of improving dynamic observables is the recent work on EDS ab initio molecular dynamics (AIMD) simulations of water[10]. DFT water with the BLYP exchange functional poorly represents water structure as seen in Figure 3 (black line). It is over-structured and has water self-diffusion coefficients that are too high (0.005-0.005 Å2 / ps)[26] compared with the experimental value of 0.23 Å2[27]. EDS was used to improve the coordination number of the oxygen-oxygen (Ow-Ow) water molecules and resulted in near perfect agreement with experimental scattering data (Figure 3). A number of unrelated observables improved as well, including RDFs and the water self-diffusion coefficient, which increased to 0.06 $\pm\,$ Å2/ps. White et al. [10] further demonstrated that the EDS bias could be transferred to excess proton-water simulations and that when combined with DFT dispersion corrections[28], the agreement further improves. The accuracy improvement with the dispersion corrections shows that EDS is general in its applicability to DFT methods.

Cortina et al. [29] used EDS to study the KPC-2 carbapenemase enzyme, which is responsible for drug resistance in the majority of carbapenem-resistant Gram-negative bacteria[30]. They used EDS to modify protein-protein distances while doing a committor analysis[31] to identify transition states in the carbapenem-enzyme complex. Cesari et al. [24] used a method similar to EDS (modified update step) to improve agreement with experimental NMR 3J coupling data for RNA oligonucleotids. After biasing, they used the simulation results to improve the underlying Amber force field and thus create a transferrable model. Cesari et al. [24] also developed a novel approach to account for experimental data uncertainty by adding auxiliary variables to the EDS update step.

EDM has been used less than EDS due to the rarity of experimental data giving an exact probability distribution. One example explored in the original EDM paper was to construct a mean field bias that mimics an alanine dipeptide being in the backbone of a protein[7]. This was done by computing a potential of mean force (PMF) for $\phi,\psi$ dihedral angles from PDB crystallography data for the alanine-alanine sequence in protein structures. This corresponds to a probability distribution and a molecular dynamics simulation of alanine dipeptide was done with its $\phi,\psi$ dihedral angles biased to match the desired values. The result was a molecular dynamics simulation where the dipeptide was biased to behave as part of a longer protein structure. A similar approach was later used by Hocky et al. [21] with EDS to built a pseudo-mean field model for actin filaments.

Another example of EDM can be found in Gil-Ley et al. [13] who used it (called targeted metadynamics) to improve agreement of RNA oligonucleatides’ dihedral angles with PDB crystallography data. They similarly built dihedral angle PMFs using the crystallography data and biased molecular dynamics simulations of RNA to adopt the dihedral angle probability distribution functions. Gil-Ley et al. [13] then took the converged EDM bias and transferred it to a different simulation of an RNA tetramer and found little to no improvement of agreement with the crystallography data. This may have been due to the underlying assumption that a PMF derived from crystallographic data is not representative of the room temperature ensemble. Nevertheless, EDM is a promising technique for incorporating experimental data as a type of auxiliary potential energy function in simulations.

3.1 Coarse-Grain Modeling

EDS and EDM are well-suited for coarse-grain (CG) modeling because they guarantee that a CG model matches experimental data. Dannenhoffer-Lafage et al. [32] showed that EDS can be applied before force-matching (a CG method) and the improvements of agreement with experimental data is maintained. Dannenhoffer-Lafage et al. [32] began with all atom molecular dynamics simulations of Ethylene carbonate with the force field from Masia et al. [33], which is known to not match center-of-mass coordination numbers from more accurate DFT calculations from Borodin and Smith [34]. Dannenhoffer-Lafage et al. [32] biased the all-atom simulations with EDS to match these known coordination numbers. They then used force-matching to create multiple CG models with one, two, and three site CG beads, using either biased all-atom or unbiased all-atom simulation data. The biased all-atom simulations showed better agreement with coordination numbers, indicating that the improvement in all-atom simulations translates to improvement in the CG model with EDS. Of course, EDS could be used directly on the CG model but that can lead to complications in how observables are calculated on CG models[35].

EDS has been applied to CG simulations of the G- and F-actin proteins as monomers and trimers[21]. Actin proteins are globular proteins that can polymerize into long semi-flexible filaments. Their hypothesis was that they could model a subsystem from within the polymerized actin structure by incorporating information about structural fluctuations from simulations of the larger system via EDS. By biasing the first and second moment of two important collective variables in actin monomer structure, they were able to observe filament-like conformations, and importantly, the fluctuations of these observables for an actin monomer, even in the system which contained only a single monomer solvated in water. Hocky et al. [21] also studied a number of questions about the EDS method in their system and found the following conclusions: (1) the linear term of EDS better matches target values and maintains system fluctuations than harmonic biases; (2) replacing the variance term in Equation 8 with a covariance matrix for all biased dimensions has faster convergence; (3) Levenberg–Marquardt[36] converges faster than stochastic gradient descent. This paper also derived a simple equation that can be used to guess the value of $\lambda$ without doing a stochastic minimization (whose accuracy depends on the distance of the target observables from the unbiased observables), and hence can serve as a good initial guess for starting an EDS simulation.

3.2 Enhanced Sampling

EDM may require enhanced sampling if there are slow degrees of freedom orthogonal to the biased collective variable. This is most conveniently treated via the extensive literature on enhanced sampling with metadynamics, since EDM is a type of metadynamics and typically implemented within a metadynamics code. EDS is not as simple because it requires that the $\left<\cdot\right>_{\tau}$ term in Equation 8 be taken over an NVT ensemble. This limits the enhanced sampling techniques to those that still give correct NVT ensemble averages. One example is parallel-tempering replica-exchange. It is possible to use metadynamics if an appropriate estimator[37] is done to compute the averages but this has not been explored in practice.

Amirkulova and White [38] demonstrated the use of enhanced sampling and EDS with the parallel-tempering well-tempered ensemble (PT-WTE)[39]. The PT-WTE method is a enhancement of enhanced parallel-tempering replica-exchange that improves exchange rates and reduces the required number of replicas[40]. PT-WTE satisfies the observable that the ensemble averages, $\left<\cdot\right>_{\tau}$ , can be computed during the course of the simulation because PT-WTE only changes the magnitude of potential energy fluctuations, not their expectations. One apparent drawback is that the method loses the one-replica observable of EDS that allows computation of dynamic observables. However, this only applies to the adaptive phase. During the fixed-bias second phase of EDS, one replica can again be used to allow analysis of dynamic observables.

Amirkulova and White [38] studied the GYG peptide with the EDS plus enhanced sampling approach. Eight simulations were conducted, each with 1, 8, or 16 replicas. The simulations had EDS, PT-WTE, and/or parallel-tempering. EDS was used to bias proton chemical shifts to improve agreement with experimental NMR data. The simulation results showed that PT-WTE improves sampling in the EDS method and does not change agreement with experimental data. PT-WTE also converged the EDS bias with fewer replicas than the PT method. One consideration in all EDS simulations is the effect on unbiased observables. The Ramachandran plots of the simulations with and without EDS bias are compared in Figure 4. When EDS and enhanced sampling are used (Figure 4 a)), the simulation explores a larger region of configurational space relative to the control simulation (Figure 4 b)). Also, the global minimum changes when using EDS in Figure 4, bringing it closer to what was found in Ting et al. [41] for the GYG sequence.

4 The Role of Pressure

Thus far we have discussed NVT and NVE ensembles. EDM and EDS add new forces to the simulation and thus affect the system virial. This can lead to undesirable density changes when the bias is subsequently used in an NPT simulation. It is possible to rectify this change in virial, and thus density, by adding a further constraint that the virial be unchanged while still maximizing entropy. As shown in the appendix, this leads to an unsolvable equation and thus it is not possible to simultaneously set the virial and other constraints. One intuitive reason for this is that the virial theorem requires the average virial be proportional to average temperature. Therefore if we tried to bias an ensemble so that its virial is fixed, then the average temperature would change, violating our constant average temperature.

One way around this challenge lies in the correlation between biased observables in EDS. When biasing multiple observables which are correlated, there is in fact no longer a unique set of Lagrange multipliers[16]. These extra degrees of freedom in the choice of Lagrange multipliers that maximize Equation 3 enable us to choose the Lagrange multipliers that minimally change the virial. Equation 8 can be modified so that our coupling constants minimize the additional pressure which they exert[42]:

[TABLE]

where $\nu$ is a parameter that controls the importance of minimizing the virial, $\Delta p_{\tau}$ is the change in the system virial pressure due to the EDS bias, $\frac{\partial s}{\partial v_{ij}}$ is the partial derivative of the observable with respect to one component of the triclinic box matrix, and $\frac{\partial v_{ij}}{\partial V}$ can be computed from the adjugate of the triclinic box matrix. When biasing multiple observables, $\Delta p_{\tau}$ should be the total over all observables. $g_{\tau}$ is in units of $s^{2}$ per energy, so $\nu$ must be in units of volume squared times $s^{2}$ per energy squared. Practically we can choose a unitless $\nu^{*}$ , which is defined from

[TABLE]

Equation 15 was used in a molecular dynamics simulation of 128 modified SPC/E water molecules. This modified SPC/E water has increased charges ( $q_{O}=-0.94$ ) to distort its coordination number. EDS was then applied to correct the coordination number using experimentally derived coordination numbers from Skinner et al. [25] with varying strengths of the virial correction. The coordination number moment definition may be found in White and Voth [6]. Figure 6 shows that coordination number and its moments are still correctly biased with the virial term. The EDS parameters were a range ( $A$ ) of 50 kj/mol, a period of 25 fs, and the Levenberg–Marquardt optimization procedure. A Nose-Hoover thermostat with a time constant of 25 fs and a timestep of 0.5 fs were used for the molecular dynamics in the LAMMPS simulation engine[43]. The NPT barostat was Parrinello-Rahman.

Figure 5 shows the effect of the virial penalty term on the EDS virial contribution. The virial contribution plotted and computed here is the mean per-particle virial energy contribution. That is $\Delta p_{\tau}/\rho/kT$ which removes the effect of particle number and energy scale. Panel a shows the virial contribution of each collective variable with $\nu^{*}=0$ (no virial minimization). The net average virial contribution is 39.2 kT. Panel b shows $\nu^{*}=10$ and there is a much lower contribution of 5.02 kT. Panel c compares the net virial contribution of three different $\nu^{*}$ values.

Figure 6 shows the impact on the convergence of the biased coordination number and its moments with the virial minimizing terms. Interestingly, adding a virial minimization term actually improves convergence. This is due to the well-known effect of regularization, which improves convergence in optimization[44]. The virial minimization term is proportional to $\alpha_{\tau}$ , so it brings down the magnitude of $\alpha_{\tau}$ . This prevents the large swings seen in the $\nu=0$ system. Thus, adding the virial minimization term not only reduces the virial contribution but can improve convergence by minimizing magnitudes of the coupling constants. The true test can be seen in the NPT portion past 20ps where the bias is fixed but the box dimensions are free. EDS with $\nu^{*}=0$ gives poor densities in NPT with a bias computed in NVT. The virial minimization reduces the change in density with a trade-off in match to the bias set-points. Virial minimization is in general recommended due to its enhanced convergence and better representation of virial forces.

5 Implementations

There are three implementations of EDS available: Colvars[45], Plumed1.3[46], and Plumed2[47]. The Colvars and Plumed1.3 implementations match the original EDS manuscript. The Plumed2 implementation is actively maintained as a plugin[48] and has features from more recent EDS articles[21]. For example, the covariance term in Equation 8 can be computed using the full sample covariance matrix or the sample variance. The update steps can also be computed using the Levenberg–Marquardt method[36, 21], instead of Equation 7. The Plumed2 implementation is recommended.

EDM is implemented in Plumed2 as well, under the “targeted metadynamics” keyword within the metadynamics biasing class. One of the common use-cases for EDM is to bias RDFs, which is not exactly a collective variable. RDFs are “function collective variables” in the sense that there is a distribution at each timestep. Thus the bias update equation, shown in Equation 10, applies to each observed point in the RDF at each step (see White et al. [7] for equations). This can lead to thousands of updates per timestep on even small systems. To improve performance and scaling, we have created an implementation in LAMMPS that is more tightly integrated with the simulation engine than Plumed2 as described in White et al. [7]. Note that the addition of hills is parallelized here, which is unusual and enables the scaling as a function of particle size. Hill addition is 90% of the CPU utilization. This was the implementation used in White et al. [7] for an ethylene carbonate electrolyte simulation and Lennard-Jones benchmark systems.

6 Conclusions

Minimal biasing techniques are an emerging area that improve the accuracy of molecular simulation and better utilize comparable experimental data. EDS and EDM are maximum entropy techniques that minimally bias an ensemble so that average observables or probability distributions of observables match set values. EDS and EDM both converge the potential to Equation 5, allowing an NVE simulation with fixed potential energy in a single replica so that dynamic observables can be computed. EDS and EDM can bias multiple observables simultaneously and be combined with enhanced sampling methods. It is possible to treat uncertainty in experimental data using methods from Cesari et al. [24]. Explicitly setting pressure is not possible in maximum entropy biasing, but in EDS it is possible to minimize the change to pressure when biasing multiple collective variables. This minimization can actually improve convergence of EDS by acting as a regularization. A variety of example systems have been presented here and there are implementations available in most simulation engines for both EDS and EDM.

Acknowledgement(s)

The authors thank Prof. Pengfei Huo for reviewing the appendix mathematics, Rainier Barrett for help preparing the manuscript and Prof. Glen Hocky for assistance in surveying the literature and preparing the manuscript.

Funding

This work was supported by the National Science Foundation CBET Div Of Chem, Bioeng, Env, & Transp Sys (Grant #1751471).

7 References

Appendix A Minimal Biasing Pressure Derivation

We would like to constrain the virial contribution to pressure while minimizing relative entropy in a biased ensemble. The constraints are:

[TABLE]

which enforces normalization. The other constraint is to set the virial pressure, $V[P^{\prime}(\vec{r})]$ , to be equal to the scalar $p_{e}$ .

[TABLE]

where $V$ is the volume. The quantity to be minimized is relative entropy of the biased probability distribution with respect to the unbiased probability distribution:

[TABLE]

To solve this system, we compute the functional derivative of the Lagrangian with respect to the biased probability distribution $P^{\prime}(\vec{r})$ . This requires a functional derivative of Equation 19, which is:

[TABLE]

Some intuition for Equation 21 can be gained from the definition of the functional derivative. Consider adding an arbitrary function $\phi(\vec{r})$ to the probability distribution scaled by $\epsilon$ , where $\epsilon$ is small. The change in the virial with this additional function is

[TABLE]

The right-hand side is independent of the probability distribution, so that changing the probability distribution decreases the virial by $N$ regardless of where the change occurs. With the functional derivative of the virial being constant, the functional derivative of the Lagrangian is

[TABLE]

which has no solution for $P^{\prime}(\vec{r})$ because $\lambda_{1}$ and $\lambda_{2}$ are colinear. Thus there is no way to constrain pressure.

Bibliography48

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Lindorff-Larsen et al. [2010] Kresten Lindorff-Larsen, Stefano Piana, Kim Palmo, Paul Maragakis, John L. Klepeis, Ron O. Dror, and David E. Shaw. Improved side-chain torsion potentials for the Amber ff 99SB protein force field. Proteins. , 78(8):1950–1958, 2010. ISSN 08873585. 10.1002/prot.22711 . · doi ↗
2Shivakumar et al. [2010] Devleena Shivakumar, Joshua Williams, Yujie Wu, Wolfgang Damm, John Shelley, and Woody Sherman. Prediction of absolute solvation free energies using molecular dynamics free energy perturbation and the opls force field. J. Chem. Theory Comput. , 6(5):1509–1519, 2010. ISSN 15499626. 10.1021/ct 900587 b . · doi ↗
3Sormanni et al. [2017] Pietro Sormanni, Damiano Piovesan, Gabriella T Heller, Massimiliano Bonomi, Predrag Kukic, Carlo Camilloni, Monika Fuxreiter, Zsuzsanna Dosztanyi, Rohit V Pappu, M Madan Babu, Sonia Longhi, Peter Tompa, A Keith Dunker, Vladimir N Uversky, Silvio C E Tosatto, and Michele Vendruscolo. Simultaneous quantification of protein order and disorder. Nat. Chem. Biol. , 13(4):339–342, 2017. ISSN 1552-4450. 10.1038/nchembio.2331 . · doi ↗
4Bonomi et al. [2017] Massimiliano Bonomi, Gabriella T. Heller, Carlo Camilloni, and Michele Vendruscolo. Principles of protein structural ensemble determination. Curr. Opin. Struct. Biol. , 42:106–116, 2017. ISSN 1879033 X. 10.1016/j.sbi.2016.12.004 . · doi ↗
5Olsson et al. [2013] Simon Olsson, Jes Frellsen, Wouter Boomsma, Kanti V Mardia, and Thomas. Hamelryck. Inference of structure ensembles of flexible biomolecules from sparse, averaged data. P Lo S One , 8(11):e 79439, 2013. ISSN 1932-6203. 10.1371/journal.pone.0079439 . · doi ↗
6White and Voth [2014] Andrew D White and Gregory A Voth. Efficient and minimal method to bias molecular simulations with experimental data. J. Chem. Theory Comput. , 10(8):3023–3030, 2014.
7White et al. [2015] Andrew White, James Dama, and Gregory a. Voth. Designing Free Energy Surfaces that Match Experimental Data with Metadynamics. J. Chem. Theory Comput. , pages 2451–2460, 2015. ISSN 1549-9618. 10.1021/acs.jctc.5b 00178 . · doi ↗
8Beauchamp et al. [2014] Kyle A Beauchamp, Vijay S Pande, and Rhiju. Das. Bayesian Energy Landscape Tilting: Towards Concordant Models of Molecular Ensembles. Biophys. J. , 106(6):1381–1390, 2014. ISSN 0006-3495. 10.1016/j.bpj.2014.02.009 . · doi ↗