Molecular simulations minimally restrained by experimental data

TL;DR

This paper derives equations for restrained molecular simulations that minimally perturb unbiased distributions while fitting experimental data, improving the integration of experimental constraints into simulations.

Contribution

It introduces a new method to restrain simulations, ensuring minimal deviation from unbiased distributions while matching experimental observable values.

Findings

Provides equations for equilibrium distributions under restraints

Suggests a method for minimally perturbing simulations

Enables accurate reproduction of experimental data within uncertainties

Abstract

One popular approach to incorporating experimental data into molecular simulations is to restrain the ensemble average of observables to their experimental values. Here I derive equations for the equilibrium distributions generated by restrained ensemble simulations and the corresponding expected values of observables. My results suggest a method to restrain simulations so that they generate distributions that are minimally perturbed from the unbiased distributions while reproducing the experimental values of the observables within their measurement uncertainties.

Equations

$$p_N(\{{\mathbf{x}}_{i=1,2,\dots ,N}\}|{\mathbf{D}}_e)$$

$$p_N(\{{\mathbf{x}}_i\};\mathbf{K},{\mathbf{D}}_e)\propto e^{-{\sum }_i\beta E({\mathbf{x}}_i)-\frac{1}{2}{(\bar{\mathbf{D}}-{\mathbf{D}}_e)}^t\cdot \mathbf{K}\cdot (\bar{\mathbf{D}}-{\mathbf{D}}_e)},$$

$$p_{\lambda }(\mathbf{x})=Z_{\lambda }^{-1}\exp (-\beta E(\mathbf{x})+{\lambda }^t\cdot \mathbf{D}(\mathbf{x}))$$

$$Z_N(\mathbf{K},\mathbf{D})$$

$$\underset{\begin{array}{c}\bar{\mathbf{D}}s.t.\\ e^{-{(\bar{\mathbf{D}}-{\mathbf{D}}_{\lambda })}^t\cdot \mathbf{K}\cdot (\bar{\mathbf{D}}-{\mathbf{D}}_{\lambda })}>ϵ\end{array}}{\max }{(\bar{\mathbf{D}}-⟨\mathbf{D}{⟩}_{\lambda })}^t\cdot N{\mathbf{\Sigma }}^{-1}\cdot (\bar{\mathbf{D}}-⟨\mathbf{D}{⟩}_{\lambda })<\alpha$$

$$Z_N(\mathbf{K},\mathbf{D})$$

$$\mathbf{K}^{\ast}=\left(\begin{array}[]{ll}\mathbf{K}&\mathbf{0}_{M\times M^{\prime}}\\ \mathbf{0}_{M^{\prime}\times M}&\mathbf{0}_{M^{\prime}\times M^{\prime}}\end{array}\right),$$

$$⟨{\mathbf{D}}^{\ast }{⟩}_{N,\mathbf{K},{\mathbf{D}}_e}$$

$$\mathbf{\Sigma}_{\lambda}^{\ast}=\left(\begin{array}[]{cc}\mathbf{\Sigma}_{\lambda}&\mathbf{c}_{0,\lambda}^{t}\\ \mathbf{c}_{0,\lambda}&\mathbf{\Sigma}_{00,\lambda}\end{array}\right)$$

$$⟨\mathbf{D}({\mathbf{x}}_1){⟩}_{N,\mathbf{K},{\mathbf{D}}_e}$$

$$⟨{\mathbf{D}}_0({\mathbf{x}}_1){⟩}_{N,\mathbf{K},{\mathbf{D}}_e}$$

$$p_N({\mathbf{x}}_1;\mathbf{K},\mathbf{D})$$

$$p_N(\mathbf{x};\mathbf{K},{\mathbf{D}}_e)$$

$$1\ll N\ll \min (\text{eigenvalues of }\mathbf{K}\cdot {\mathbf{\Sigma }}_{\lambda }).$$

$$p_N({\mathbf{x}}_1;\mathbf{K},{\mathbf{D}}_e)\propto p_0({\mathbf{x}}_1)e^{(\lambda +{\mathbf{\Sigma }}_{\lambda }^{-1}\cdot ({\mathbf{D}}_e-⟨\mathbf{D}{⟩}_{\lambda }){)}^t\cdot (\mathbf{D}({\mathbf{x}}_1)-{\mathbf{D}}_e)},$$

$${\Delta }_{\text{KL}}[p_{{\lambda }^{\prime }}||p_N]$$

$$e^{G(\mathbf{D}(\mathbf{x}),N)}$$

$$p(x,y;f_a,f_b)=e^{-E(x,y;f_a,f_b)}$$

$${\lambda }_j={\lambda }_{j-1}+{\mathbf{\Sigma }}_{{\lambda }_{j-1}}^{-1}\cdot ({\mathbf{D}}_e-⟨\mathbf{D}{⟩}_{{\lambda }_{j-1}})$$

Full text

I Introduction

Molecular simulations are routinely used to generate models of molecular structures and their equilibrium distributions Dror et al. (2012). They provide structural and dynamic information at atomistic resolution that are critical to the interpretation of experimental data Bottaro and Lindorff-Larsen (2018). Inaccuracy in the empirical force fields used to compute the potential energy of the biomolecules in the simulations, however, can cause the simulations to generate inaccurate equilibrium distributions of molecular structures, manifested by the disparity between the experimentally measured values of observables and the predictions by the simulations Gunsteren, Dolenc, and Mark (2008). Biasing the simulations toward reproducing the experimental values of the observables is a promising approach to improving the accuracy of the simulations Boomsma et al. (2014); Perez et al. (2016). How to introduce minimal bias, however, remains an open question, and no such method exists that takes into account the measurement uncertainties associated with the experimental data.

Because many experimental techniques, such as nuclear magnetic resonance (NMR) and small-angle X-ray scattering (SAXS), measure the ensemble average of observables, a natural and popular method to incorporate such experimental data is to simultaneously simulate many replicas of the molecular system while restraining the average of the instantaneous values of any observable of interest to be close to its experimental value Lindorff-Larsen et al. (2005); Vendruscolo (2007); Roux and Islam (2013); Hummer and Köfinger (2015). I will refer to such replica-based restrained simulations simply as restrained ensemble simulations. The equilibrium distribution generated by restrained ensemble simulations depends on the number of replicas and the strength of the restraints, the latter of which should reflect the measurement uncertainties of the observables Hummer and Köfinger (2015). There lacks, however, a clear theoretical underpinning for the choice of the number of replicas Olsson and Cavalli (2015).

The restrained simulations should be no more biased by the restraints than is necessary, which is a principle that motivated the alternative maximum entropy method, which seeks a distribution that reproduces the experimental values of the observables and is minimally perturbed from the unbiased distribution  Boomsma, Ferkinghoff-Borg, and Lindorff-Larsen (2014); Różycki, Kim, and Hummer (2011); Groth et al. (1999); Massad et al. (2007). It has been established that in some limit of infinite number of replicas and infinitely strong restraints, restrained ensemble simulations generate a statistically equivalent distribution as the maximum entropy method Pitera and Chodera (2012); Roux and Weare (2013); Cavalli, Camilloni, and Vendruscolo (2013); Olsson and Cavalli (2015). In practice, however, the restraining strength is finite, in which case the two methods are not equivalent. In contrast to restrained ensemble simulations, the maximum entropy method constrains every observable to its exact experimental value, thus it does not take into account the measurement uncertainties, which are ubiquitous in all experiments.

Here I present theoretical results that enable the determination of the largest number of replicas to use in restrained ensemble simulations with fixed strength of restraints so that they generate a distribution minimally perturbed from the unbiased distribution while reproducing the experimental values of the observables within their measurement uncertainties. As I will demonstrate below, this is equivalent to the determination of the weakest strength of restraints to use at fixed number of replicas, provided that the latter is sufficiently large. The key result is an equation (Eq. LABEL:eqn:D-expected) that predicts the expected value of any–restrained and unrestrained–observable in restrained ensemble simulations for different numbers of replicas, thus permitting the selection of the appropriate number without computationally expensive trial and error. In addition, I derive an equation (Eq. 12) for the distribution generated by restrained ensemble simulations with a finite number of replicas and a finite restraining strength, which can be used to perform reweighted simulations that take into account of the experimental data with their attendant measurement uncertainties. The validity of these theoretical results are demonstrated numerically by Monte Carlo simulations. I illustrate the application of these results using the molecular dynamics (MD) simulations of the polyalanine peptide Ala5 restrained by experimental nuclear magnetic resonance (NMR) data.

II Theory

Consider an ensemble consisting of $N$ replicas of the molecular system, each with the atomic coordinates $\mathbf{x}_{i=1,2,\dots,N}$. Experiments measure the ensemble average of the individual contributions of the observables: $\bar{\mathbf{D}}=N^{-1}\sum_{i}\mathbf{D}_{i}$, where $\mathbf{D}_{i}\equiv\mathbf{D}(\mathbf{x}_{i})=(D_{1}(\mathbf{x}_{i}),D_{2}(\mathbf{x}_{i}),\dots,D_{M}(\mathbf{x}_{i})^{t}$ are the instantaneous values of the vector of $M$ observables for the $i$’th replica. In an unbiased simulation, $\bar{\mathbf{D}}$ may deviate from experimental values $\mathbf{D}_{e}$. The posterior distribution given the experimental data Rieping, Habeck, and Nilges (2005); Hummer and Köfinger (2015) is

[TABLE]

where $p_{0}(\mathbf{x}_{i})=Z_{0}^{-1}\exp(-\beta E(\mathbf{x}_{i}))$ is the Boltzmann distribution of the molecular configurations under the potential energy function $E$. Assuming a normal distribution in each experimental value, $p(\mathbf{D}_{e}|\{\mathbf{x}_{i}\})=(\prod_{\alpha}2\pi\delta_{\alpha}^{2})^{-M/2}\exp(-\frac{1}{2}{(\bar{\mathbf{D}}-\mathbf{D}_{e})}^{t}\cdot{\mathbf{K}}\cdot{(\bar{\mathbf{D}}-\mathbf{D}_{e})})$, where $\mathbf{K}=\text{diag}(\delta_{1}^{-2},\delta_{2}^{-2},\dots,\delta_{M}^{-2})$ is the inverse of the statistical variances of the experimental measurements, $p_{N}(\{\mathbf{x}_{i}\}|\mathbf{D}_{e})$ becomes

[TABLE]

which is the distribution in the restrained ensemble simulations. The strength of the restraints, $\mathbf{K}$, is thus related to the measurement uncertainties associated with the experimental data. In the restrained ensemble simulations, $N$ replicas of the molecular system are simulated with the same energy function and at the same temperature, and they interact with one another only through the quadratic term in the exponent in Eq. 2.

A similar equation was derived by Hummer and Köfinger (Eq. 26 of Ref. Hummer and Köfinger, 2015), with the difference that they scaled the restraining term by $N$ and introduced a parameter ($\theta$) expressing the level of confidence in the unbiased reference distribution. My results below show that these two equations are equivalent at large $N$.

For a fixed $\mathbf{K}$, the marginal distribution for a single replica, $p_{N}(\mathbf{x}_{1};\mathbf{K},\mathbf{D})=\int\prod_{i=2}^{N}p_{N}(\{\mathbf{x}_{i}\};\mathbf{K},\mathbf{D}_{e})$, tends to the unbiased distribution $p_{0}(\mathbf{x}_{1})$ with increasing $N$. Larger $N$ thus corresponds to higher confidence in the unbiased distribution. A reasonable choice of $N$ is to take the largest value such that the expected values of the observables in the restrained ensemble simulations remain within the measurement uncertainties of the experimental values: this generates a distribution minimally perturbed from the unbiased distribution while still reproducing the experimental data.

I will analyze the behavior of restrained ensemble simulations at sufficiently large numbers of replicas, i.e. $N\gg 1$. Let

[TABLE]

be a biased distribution for an arbitrary vector $\lambda$, which becomes the distribution generated by the maximum entropy method if $\lambda$ is the corresponding Lagrangian multiplier Roux and Weare (2013). Dropping the subscript from the constant vector $\mathbf{D}_{e}$, the partition function for the restrained ensemble of $N$ replicas is

[TABLE]

where $\langle\cdots\rangle_{N,\lambda}$ signifies the ensemble average according to the probability distribution of $N$ independent replicas (i.e., $\mathbf{K}=\mathbf{0}$), each sampled according to the biased distribution $p_{\lambda}(\mathbf{x})$, and $\rho_{N,\lambda}(\bar{\mathbf{D}})$ is the corresponding probability distribution of the random variable $\bar{\mathbf{D}}=N^{-1}\sum_{i}\mathbf{D}_{i}$. In this ensemble of $N$ independent replicas, $\{\mathbf{D}_{i}\}$ are independent and identically distributed (i.i.d.) random variables. Acccording to the central limit theorem, the random variable $\bar{\mathbf{D}}=N^{-1}\sum_{i}\mathbf{D}_{i}$, for sufficiently large $N$, tends to the normal distribution $\rho_{N,\lambda}(\bar{\mathbf{D}})\approx(2^{M}\pi^{M}|\mathbf{\Sigma}_{\lambda}|/N)^{-1/2}e^{-\frac{1}{2}{(\bar{\mathbf{D}}-\langle\mathbf{D}\rangle_{\lambda})}^{t}\cdot{N\mathbf{\Sigma}^{-1}}\cdot{(\bar{\mathbf{D}}-\langle\mathbf{D}\rangle_{\lambda})}}$, where $\langle\mathbf{D}\rangle_{\lambda}=\int d\mathbf{x}\,\mathbf{D}(\mathbf{x})p_{\lambda}(\mathbf{x})$ are the expected values of the observables $\mathbf{D}(\mathbf{x})$ in the biased distribution (I will use $\langle\cdots\rangle_{\lambda}$ to denote the expected values in the biased distribution with parameter $\lambda$), $\mathbf{\Sigma}_{\lambda}=\langle(\mathbf{D}(\mathbf{x})-\langle\mathbf{D}\rangle)\cdot(\mathbf{D}(\mathbf{x})-\langle\mathbf{D}\rangle)^{t}\rangle_{\lambda}$ is the corresponding covariance matrix, and $|\mathbf{\Sigma}_{\lambda}|$ is the determinant of matrix $\mathbf{\Sigma}_{\lambda}$. If $\mathbf{D}_{\lambda}\equiv\mathbf{D}-N\mathbf{K}^{-1}\cdot\lambda$ and $\langle\mathbf{D}\rangle_{\lambda}$ are close, in the sense that

[TABLE]

for a sufficiently small $\epsilon>0$ and some value of $\alpha\geq 3$, such that the normal distribution is a good approximation to $\rho_{N,\lambda}(\bar{\mathbf{D}})$ for all $\bar{\mathbf{D}}$ around $\mathbf{D}_{\lambda}$, $\bar{\mathbf{D}}$ can be integrated out, and the partition function in Eq. 4 becomes

[TABLE]

where $I$ is the identity matrix. Eq. 6 forms the basis for deriving the key results of this work.

My first key result is that the expected value of any observable, restrained or unrestrained, in the restrained ensemble simulations can be predicted from the biased simulation of $p_{\lambda}(\mathbf{x})$. Let $\mathbf{D}^{\ast}=(D_{1},D_{2},\dots,D_{M},D_{M+1},\dots,D_{M+M^{\prime}})^{t}$ be the vector of restrained and unrestrained observables, where $\mathbf{D}=(D_{1},D_{2},\dots,D_{M})$ are the $M$ restrained observables and $\mathbf{D}_{0}=(D_{M+1},\dots,D_{M+M^{\prime}})^{t}$ are the $M^{\prime}$ unrestrained observables. Introducing

[TABLE]

$\mathbf{D}_{e}^{\ast}=(\mathbf{D}_{e}^{t}\ \,\mathbf{0}_{1\times M^{\prime}})^{t}$, and $\lambda^{\ast}=(\lambda^{t}\ \,\mathbf{0}_{1\times M^{\prime}})^{t}$, the expected values of the observables $\mathbf{D}^{\ast}(\mathbf{x})$ in the restrained ensemble simulations are given by

[TABLE]

where $\mathbf{\Sigma}^{\ast}_{\lambda}=\langle(\mathbf{D}^{\ast}-\langle\mathbf{D}^{\ast}\rangle)\cdot(\mathbf{D}^{\ast}-\langle\mathbf{D}^{\ast}\rangle)^{t}\rangle_{\lambda}$ is the covariance between all the observables in the biased simulation.

Writing $\mathbf{\Sigma}^{\ast}$ in the block form

[TABLE]

where $\mathbf{\Sigma}_{00,\lambda}$ is the $M^{\prime}\times M^{\prime}$ covariance matrix for the unrestrained observables $\mathbf{D}_{0}$ and $\mathbf{c}_{0,\lambda}$ is the $M^{\prime}\times M$ covariance matrix between the unrestrained and the restrained observables, the expected values of the restrained observables and the unrestrained observables are, respectively,

[TABLE]

and

[TABLE]

Eq. LABEL:eqn:restrained-D-expected enables the determination of the largest $N$ such that the expected values of the restrained observables remain within the measurement uncertainties of their experimental values.

My second key result is that the marginal distribution for a single replica generated by the restrained ensemble simulations is

[TABLE]

where I used Eq. 6 for both $Z_{N}(\mathbf{K},\mathbf{D})$ and $Z_{N-1}(\frac{(N-1)^{2}}{N^{2}}\mathbf{K},\frac{N\mathbf{D}-\mathbf{D}_{1}}{N-1})$, and $N\approx N-1$ for large $N$.

For $N\gg 1$, $\frac{1}{2N}{(\mathbf{D}_{1}-\mathbf{D})}^{t}\cdot{(N\mathbf{K}^{-1}+\mathbf{\Sigma}_{\lambda})^{-1}}\cdot{(\mathbf{D}_{1}-\mathbf{D})}\approx 0$, and Eq. 12 reduces to

[TABLE]

which depends on $N$ through the product $N\mathbf{K}^{-1}$. This suggests that a restrained ensemble simulation with $N$ replicas and the strength of restraints $\mathbf{K}$ is equivalent to another one with $N^{\prime}$ replicas and the strength of restraints $\mathbf{K}^{\prime}=(N^{\prime}/N)\mathbf{K}$. This equivalency, to my knowledge, has not been previously demonstrated except for the harmonic potentials Roux and Weare (2013), and it mathematically justifies the scaling by $N$ of the restraining term proposed by Hummer and Köfinger Hummer and Köfinger (2015).

Eq. 12 suggests that the restrained ensemble simulations generate the same distribution as the maximum entropy method if the number of replicas satisfies the following inequality

[TABLE]

This condition is the generalization of the condition $N\ll K$ derived for the special case of harmonic potentials Roux and Weare (2013). Under this condition, $(I+N\mathbf{K}^{-1}\mathbf{\Sigma}_{\lambda}^{-1})^{-1}\approx I$, the corresponding single-replica distribution becomes

[TABLE]

which is the same as $p_{\lambda}(\mathbf{x}_{1})$ if $\lambda$ takes the value of the Lagrangian multiplier $\lambda_{\text{ME}}$ in the maximum entropy method such that $\langle\mathbf{D}\rangle_{\lambda_{\text{ME}}}=\mathbf{D}_{e}$. Under this condition, Eq. LABEL:eqn:restrained-D-expected becomes $\langle\mathbf{D}(\mathbf{x}_{1})\rangle_{N,\mathbf{K},\mathbf{D}_{e}}\approx\mathbf{D}_{e}$.

The Kullback-Leibler divergence from the distribution in the restrained ensemble simulations (Eq. 12) to that of $p_{\lambda^{\prime}}(\mathbf{x})$ is

[TABLE]

where

[TABLE]

Eq. 16 enables the quantitative estimation of the difference between the distributions generated by the restrained ensemble simulations and by the maximum entropy method. Minimization of $\Delta_{\text{KL}}[p_{\lambda_{\text{ME}}}||p_{N}]$ with respect to $N$ yields the number of replicas to use such that the restrained ensemble simulations generate a distribution closest to that of the maximum entropy method.

III Results and Discussions

I demonstrate the validity of Eq. LABEL:eqn:D-expected and Eq. 16 using the Monte Carlo simulations of a 2-dimensional potential that generates a distribution of mixed Gaussians:

[TABLE]

where $f_{a}$ and $f_{b}$ are the mixing coefficients constrained by $f_{a}+f_{b}=1$. The parameters used in the simulations are $\sigma_{a}=((1.5,0)^{t},(0,1.5)^{t})$, $(x_{a},y_{a})=(-1,0)$, and $\sigma_{b}=((1.5,0.1)^{t},(0.1,1.5)^{t})$, $(x_{b},y_{b})=(0.5,2.5)$. The restrained ensemble simulations (Eq. 2) and the maximum entropy method are applied to a potential with $f_{a}=0.3$, restraining two observables, $x$ and $\sqrt{x^{2}+y^{2}}$, to reference values that are equal to the observables’ mean values in a reference distribution that has $f_{a}=0.7$. The restraining strength is $\mathbf{K}=\text{diag}(100,100)$.

For sufficiently large $N$, Eq. LABEL:eqn:D-expected predicts the expected values of both the restrained and unrestrained observables in quantitative agreement with their mean values in the restrained ensemble simulations (Fig. 1). As $N$ increases, the distribution generated by the restrained ensemble simulations tends to the unbiased distribution ($\lim_{N\rightarrow\infty}\Delta_{\text{KL}}[p_{0}||p_{N}]=0$). Notably, $\Delta_{\text{KL}}[p_{0}||p_{N}]$ is relatively insensitive to $N$ around the largest value of $N$ for which the expected values of the restrained observables remain close to their reference values .

My results suggest the following method of setting $N$ so that the simulations are minimally restrained by experimental data: Iteratively perform a series of biased simulations of $p_{\lambda_{j}}(\mathbf{x})$, updating $\lambda_{j}$ by either

[TABLE]

or other update schemes Stinis (2005) such that $\lim_{j\rightarrow\infty}\lambda_{j}=\lambda_{\text{ME}}$. For each iteration, compute $\langle\mathbf{D}(\mathbf{x}_{1})\rangle_{N,\mathbf{K},\mathbf{D}_{e}}$ using Eq. LABEL:eqn:restrained-D-expected and select the largest $N_{j}$ so that the expected values of all the restrained observables are within the measurement uncertainties from their experimental values. Continue until $N_{j}$ no longer changes.

I illustrate this method of setting $N$ by MD simulations of the polyalanine peptide Ala5 restrained by the ${}^{3}J$ couplings measured by NMR experiments Graf et al. (2007); Wickstrom, Okur, and Simmerling (2009). All MD simulations were performed using the OpenMM software Eastman et al. (2013), with Amber99SB-ILDN force field Lindorff-Larsen et al. (2010) and a generalized Born implcit solvent model Onufriev, Bashford, and Case (2004). Langevin integrator Leimkuhler and Matthews (2015) was used to maintain the temperature of the molecular system at 300 K. When Ala5 was simulated without any restraints, there was substantial discrepancy between the predictions by the unbiased simulation and the experimental values for the ${}^{3}J$ couplings between H${}^{\text{N}}$ and Hα (denoted as ${}^{3}J_{\text{H}^{\text{N}},\text{H}^{\alpha},r}$ for $r=2,3,4,5$), which have been used to characterize the four backbone torsions $\phi_{r=2,3,4,5}$ of Ala5 Vögeli et al. (2007) (Fig. 2). Below, ${}^{3}J_{\text{H}^{\text{N}},\text{H}^{\alpha},r=2,3,4,5}$ were used as the restrained observables.

The same molecular system was then simulated with a series of biasing potentials (Eq. 3), in which $\lambda_{j}$ in the $j$’th simulation was determined by Eq. 19. Eq. LABEL:eqn:restrained-D-expected was used to predict the expected values of the restrained observables (${}^{3}J_{\text{H}^{\text{N}},\text{H}^{\alpha},r=2,3,4,5}$) in the restrained ensemble simulations for different values of $N$, using the values of $\lambda_{j}$ and the corresponding $\langle\mathbf{D}\rangle_{\lambda_{j}}$ and $\mathbf{\Sigma}_{\lambda_{j}}$ in different biased simulations as inputs. The largest $N$ ($N_{\text{max}}$) for which the predicted values of all the restrained observables are within the statistical uncertainty from the experimental values is determined for each $j$. In this case, only 1 iteration of biased simulation was necessary to have a converged estimate of $N_{\text{max}}$ (Fig. 2).

In practice, it is easier to select the number of replicas $N$ based on the available computing resources (e.g., the number of CPU cores). Once $N_{\text{max}}$ has been determined by the above procedure, one can perform the restrained ensemble simulations using $N$ replicas and the scaled strength of restraints $\mathbf{K}_{\text{min}}=(N/N_{\text{max}})\mathbf{K}$, because the restrained ensemble simulations with the parameters $(N,\mathbf{K}_{\text{min}})$ and with the parameters $(N_{\text{max}},\mathbf{K})$ generate equivalent marginal distributions (Eq. LABEL:eqn:marginal-func-of-NiK).

My results suggest that the restrained ensemble simulations at large $N$ are equivalent to a single-replica simulation of the reweighted distribution in Eq. 12 (which becomes Eq. LABEL:eqn:marginal-func-of-NiK as $N\rightarrow\infty$). In this case, $N$ can be regarded as a parameter that tunes the strength of the restraints, which are proportional to the inverse of the covariance of the experimental measurements. As $N$ tends to infinity, the resulting distribution becomes unbiased. To demonstrate this reweighting approach, the Ala5 system was simulated according to the marginal distribution in Eq. 12 for different values of $N$, using different values of $\lambda_{j}$ and the corresponding $\langle\mathbf{D}\rangle_{\lambda_{j}}$ and $\mathbf{\Sigma}_{\lambda_{j}}$ as inputs. For $j\geq 1$, the expected values of the restrained (${}^{3}J_{\text{H}^{\text{N}},\text{H}^{\alpha},r=2,3,4,5}$) and unrestrained (${}^{3}J_{\text{H}^{\text{N}},\text{C}^{\beta},r=2}$ and the radius of gyration $R_{g}$) observables in these reweighted simulations are in good agreement with the corresponding values predicted by Eq. LABEL:eqn:D-expected (Fig. 2). Such reweighted simulations may offer a simple alternative to replica-based restrained ensemble simulations as a way to incorporate into MD simulations experimental data with their associated measurement uncertainties.

IV Acknowledgments

H. X. thanks Woody Sherman for helpful suggestions on the manuscript.

V Notes

The author declares no competing financial interest.