Restrained-Ensemble Molecular Dynamics Simulations Based on Distance Histograms from Double Electron-Electron Resonance Spectroscopy
DEER (double electron electron resonance) spectroscopy is a powerful pulsed ESR (electron spin resonance) technique allowing the determination of spin-spin distance histograms between site-directed nitroxide label sites on a protein in their native environment. However, incorporating ESR/DEER data in structural refinement is challenging because the information from the large number of distance histograms is complex and highly coupled. Here, a novel restrained-ensemble molecular dynamics simulation method is developed to incorporate the information from multiple ESR/DEER distance histograms simultaneously. Illustrative tests on 3 coupled spin-labels inserted in T4 lysozyme show that the method efficiently imposes the experimental distance distribution in this system. Different rotameric states of the \xcf\x871 and \xcf\x872 dihedrals in the spin-labels are also explored by restrained ensemble simulations. Using this method, it is hoped that experimental restraints from ESR/DEER experiments can be used to refine structural properties of biological systems.
I. Introduction<!>II. Theoretical Developments<!>IV. Computational Details<!>IV. Results and Discussion<!>V. Summary
<p>DEER (Double Electron Electron Resonance) is a powerful pulsed double-quantum coherence (DQC) electron spin resonance (ESR) spectroscopy that allows the determination of the nitroxide-nitroxide distance histogram between two methanethiosulfonate spin-labels (MTSSL) attached to a protein via cysteine residues.1–4 ESR/DEER, in conjunction with site-directed spin labeling (SDSL) techniques,5,6 now allows one to measure the distance histograms between an unprecedented large number of sites, which requires new computational strategies for structural refinement. For example, up to 51 pairs of experimental ESR/DEER distance histograms have been measured for T4 lysozyme, a small benchmark protein of 165 residues (H. Mchrouab, personal communication). Exploiting the data from ESR/DEER in structural refinement, however, is confronted by several challenges. The measured histograms are information-rich and often have complex structures with multiple peaks. Furthermore, the nitroxide spin-labels (Figure 1) display considerable internal flexibility, 7 which implies that multiple rotameric states must be taken into account and that it is unrealistic to consider a single rotameric state. This uncertainty increases the difficulty in translating the measured distance histograms between spin-labels into a structural refinement restraint. Naively imposing a simple average distance between Cα backbone atoms is definitely not a valid approximation. Of particular importance, a given spin-label can be involved in multiple distance histograms, which implies that the large body of data is highly coupled. For instance, this precludes an approach where the spin-label pairs could be treated one-by-one. Furthermore, the characterization of the accessible rotamers of the spin-label side chain is difficult.7–13 In summary, the body of data from ESR/DEER distance histograms is complex, information-rich, highly coupled, and voluminous. While there have been previous efforts to utilize the information from ESR/DEER for protein structure determination and progress has been made using rotamer libraries,2–4 a general strategy is still lacking.</p><p>In familiar energy-based structural refinement algorithms, whether it is based on energy-minimization or molecular dynamics (MD) simulations, one typically introduces a biasing potential that depends on the configuration x of the system in order to enforce that a given property q(x) will match the experimental value Q. One of the simplest form is a quadratic potential of the type, k(q(x)−Q)2, with a restraint force constant k. This is the approach that is generally adopted in structural refinement based on X-ray crystallography or in nuclear magnetic resonance (NMR).14 However, by doing so we unavoidably impose also the magnitude of the fluctuations of the property q, indirectly corrupting the atomic configurations in a way that is arbitrary and undesirable. An alternative strategy is based on the notion of "restrained ensemble" MD simulations introduced by Vendruscolo and coworkers, 15–17 and exploited by Im and co-workers to determine the structure of membrane-bound peptides on the basis of solid state NMR observables.18–21 The restrained-ensemble MD simulation scheme consists in carrying out parallel MD simulations of N replicas of the basic system in the presence of a biasing potential that enforces the ensemble-average of a given property toward its known experimental value. Working from the perspective of an ensemble is appealing since it attempts to mimic the conditions under which bulk measurements are carried out experimentally. Intuitively, it is hoped that if the ensemble comprises a fairly large number of replicas, the biasing potential enforcing the ensemble-average property constitutes only a mild perturbation acting on any individual replica without causing large unwanted distortions. The idea has been illustrated by showing that solid state NMR properties within the restrained-ensemble MD rapidly converge toward unique statistical distributions as N becomes large.18–20 Recently, it was demonstrated that the restrained-ensemble MD simulation scheme is formally related to the maximum entropy method for biasing thermodynamic ensembles.22–24 This helped clarify the general significance of the restrained-ensemble MD simulations, and the underlying conditions for its validity.</p><p>Our goal with this article is to exploit these recent developments to formulate a computational approach to utilize the information from multiple ESR/DEER distance histograms in protein structure refinement. Here, an approach based on a mean-field version of the restrained-ensemble MD is elaborated and illustrated to the case of three nitroxide spin-labels attached to the small soluble protein T4 lysozyme. The article is concluded by outlining future applications of this approach in structural refinement.</p><!><p>An issue of general interest in computer simulations is to incorporate information from experiments into a structural model. An important caveat in pursuing this goal is to avoid corrupting the resulting model with spurious and arbitrary biases. Let us consider a system that is described by the potential energy U(X), where X ≡ {r1,
r2, …,
rN} represents all atomic coordinates. The equilibrium Boltzmann distribution is,</p><p>
(1)P0(X)=e-βU(X)∫dXe-βU(X) where β = 1/kBT. Without any additional biases, this model based on potential energy U(X) represents all our prior knowledge about the system. Unfortunately, this model is imperfect and, for example, the average over some property q(X),</p><p>
(2)〈q〉=∫dXq(X)P0(X) differs from the value Q known from experiment. In this regard, it is natural to attempt to incorporate the experimental information into the model. In practice, this problem can be formulated within the maximum entropy method,22 whereby one seeks to maximize an excess cross-entropy functional under the constraint that the resulting average must reproduce the experimental data.23,24 This yields a solution of the general form</p><p>
(3)P(X)=e-β[U(X)-kBTλq(X)]∫dXe-β[U(X)-kBTλq(X)] where λ is a Lagrange multiplier that must be adjusted to satisfy Eq. (2). It is as if the system is evolving on an effective potential energy surface, Ueff(X) = U(X) − kBTλq(X).</p><p>An interesting special case occurs when the set of properties corresponds to a continuous distribution along some variable ξ that is a function of X. This is the type of data that is provided by ESR/DEER spectroscopy.1 For example, ξ could be the distance between two atoms or the position of an atom along some axis z. Dividing the values of ξ into a set of L discrete bins {ξ(i)} of width Δξ, the probability in each bin i can be defined as the average of the population operator,</p><p>By construction, the set of function hi is complete and satisfy the normalization condition, Σi hi = 1. In the limit that Δξ is very small, the quantity hi(ξ)/Δξ turns into a delta function. We seek to modify the equilibrium distribution P0(X) such that the averages 〈hi〉 will match the set of experimental values Hi, for i = 1, …, L (also normalized, Σi Hi = 1). This problem can be easily formulated using the maximum cross-entropy method, which yields,</p><p>
(5)P(X)=e-β[U(X)-kBT∑iλihi(ξ(X))]∫dXe-β[U(X)-kBT∑iλihi(ξ(X))] where λi are Lagrange multipliers to be adjusted in order to match the set Hi. The system is evolving on the effective potential energy surface, Ueff = U − kBT Σi λihi. In the limit where Δξ become very small, the sum of λihi(ξ) support a continuous biasing potential Ubias as a function of the coordinate ξ,</p><p>
(6)Ubias(ξ)=-ln[H(ξ)h(ξ)] where H(ξ) is the desired experimental distribution, and h(ξ) is the distribution from the unbiased simulation. Ubias(ξ) is related to the potential of mean force (PMF), W(ξ), of the system along the coordinate ξ, with W(ξ) − W(ξ*) = −kBT ln[h(ξ)/h(ξ*)], where ξ* is an arbitrary reference point. This is true because the functions hi(ξ) form a complete and non-overlapping (orthogonal) basis set along the axis ξ. This shows that the most straightforward route to bias a system toward a single experimental distribution function consists in determining the PMF as a function of ξ, and then biasing the system using the potential Ubias defined by Eq. (6). This form is valid for any specific variable that does not involve indirect or hidden correlations, such as any internal degrees of freedom internal in a given residue, or a distance between residue i and j in a protein. This form does not, however, apply to cases where multiple degrees of freedom are superimposed to construct a single average distribution, such as the particle-particle PMF in a liquid for example.</p><p>An important extension of the above discussion arises when trying to use the information from multiple experimental ESR/DEER distance distributions to bias the system. Direct application of the maximum entropy principle becomes essentially impractical because it requires the determination of a large number of unknown Lagrange multipliers by solving iteratively until convergence is reached. An alternative strategy is the restrained-ensemble MD simulation scheme.15–21 It consists in carrying out parallel MD simulations of N replicas of the basic system in the presence of a biasing potential that restrains the ensemble-averaged property toward the experimental value Q. The restrained-ensemble MD simulation is formally equivalent to the maximum entropy principle for biasing an average property toward an experimental value.22–24 In the restrained-ensemble MD simulations, we consider an ensemble of N parallel systems s, each with the coordinates Xs and the potential energy U(Xs), and the ensemble-average of system</p><p>
(7)q¯(X1,…,XN)=1N∑s=1Nq(Xs) is restrained to match the set of experimental distributions by the potential URE</p><p>
(8)URE(X1,…,XN)=12K(q¯(X1,…,XN)-Q)2 where K is a large force constant. In the case of ESR/DEER, the experimental data is a set of distance histogram bins, Hij(r), between ij pairs of spin-labels, where r is the distance between the oxygen atoms of the nitroxide spin-labels, which carries the unpaired electrons. However, a rigorous application of the restrained-ensemble would require extracting a single count for each copy of the system, that is, the distance |
ris-rjs| between spin-labels i and j for the s-th copy. A fair amount of useful information is wasted with this strict application of the restrained-ensemble. Since ESR/DEER experiments are typically carried out for spin-labels that are well separated, by a distance on the order of 10–50 Å, it can be safely assumed that the spin-labels are uncorrelated from one another. Considering the joint probability of the s-th copy of two nitroxides i and j, 
P(ris,rjs), the ensemble average distance histogram of the ij pair is,</p><p>
(9)h¯ij(r)=1N∑s=1N∫drisdrjsδ(∣ris-rjs∣-r)P(ris,rjs), which is replaced by the mean-field approximation,</p><p>The underlying assumption supporting the ansatz 
P(ris,rjs)≈P(ris)P(rjs′) is that the conformational distribution of the s-th copy of the spin-label i is not affected by the conformational distribution of the s′-th copy of the spin-label j, which is reasonable if the two spin-labels are not in direct contact. Accordingly, the N copies of the two spin-labels yield N × N separate entries to construct the distance histogram of the ij pair. In other words, it is assumed that the distribution of conformational rotamers of spin-label i is not affected by the conformational rotamers of spin-label j. This is a reasonable approximation as long as the spin-labels i and j are well separated from one another and do not directly collide. Under this condition, the population of accessible rotameric states for spin-labels i and j are statistically independent. Such a treatment of the multiple copies corresponds to a mean-field approximation.25</p><p>In principle, the discrete counting function based on Eq. (4) could be used to construct the histograms in Eq. (10), although this leads to cumbersome discontinuities. To avoid these problems, a smooth differentiable Gaussian of width σ is used in the construction of the histograms, discretizing the distance r via r = nΔr, where Δr is the bin width. It follows that the n-th bin of the histogram for the ij pair is calculated as,</p><p>
(11)h¯ij(n)=1N2∑s=1N∑s′=1N12πσ2e-(nΔr-∣ris-rjs′∣)22σ2 where the set of copies for the first spin-label is given by s ∈ {i} and the set of copies for the second spin-label is given by ′ ∈ {j}. All possible pairs of distances between the copies of two spin-labels are collected into a single histogram, according to a mean-field approximation. The histograms, both from experiment and the restrained-ensemble simulations, are properly normalized, namely, Σ{bin
n}
h̄ij(n) Δr =Σ{bin
n}
Hij(n)Δr = 1, for all pairs ij (the Jacobian is implicitly incorporated into the integrand and the histograms have dimension of inverse length). This form has the advantage of allowing continuous first derivatives for the restraining potential. Furthermore, it is also consistent with the notion that the unpaired electron is spatially delocalized over the oxygen of the nitroxide spin-label. The restraining energy is given by,</p><p>
(12)URE=12K∑{pairij}∑{binn}(h¯ij(n)-Hij(n))2 where the force constant K has dimension of energy per length square. The first derivative of the restrained-ensemble energy with respect to the position s-th copy of the i-th spin-label, 
ris, is,</p><p>All the computational details for testing the restrained-ensemble of ESR/DEER histograms are given in the following section.</p><!><p>All the simulations were carried out using the program CHARMM.26 The atomic structure of T4 lysozyme (T4L) determined by X-ray crystallography27 was taken from the Protein Data Bank (PDB id 2LZM). T4L is a small protein of 165 residues that has often served as a model system to benchmark electron spin resonance (ESR) spectroscopy studies in conjunction with site-directed spin-labeling (SDSL) techniques. 5,6 Three residues, Glu62, Thr109, and Ala134 were mutated to cysteine and nitroxide spin-labels were attached at these positions. A novel energy restraint based on Eq. (12) was implemented in the MMFP module of CHARMM to couple multiple copies and enforce the experimental distributions from ESR/DEER. The analytical first derivative of the new energy restraint was validated by comparison with the finite-difference derivative using the TEST FIRST command. The restrained-ensemble force constant K was set to 10000 (kcal/mol)/Å2, and the natural spread σ of the Gaussian was set to 1.7 Å. This value is justified because the width of the ESR/DEER histograms is on the order of about 10 Å. The bin width of the histogram was set to 1.0 Å and the histograms run from 1 to 60 Å. Each of the spin-labels was expanded into N copies (N = 1, 10, and 25) using the multiple-copy Locally Enhanced Sampling (LES) method,25,28,29 while the remainder of the system was simulated with a single copy. The simulations with LES were setup using the REPLICA and BLOCK commands of the program CHARMM.26 In the LES simulations, the interactions of the replicated spin-labels with the rest of the system were scaled by 1/N while the internal interactions of the replicated spin-labels were unscaled. The different replicas of a given spin-label do not interact with one another. The CHARMM force field PARAM2730 with the CMAP correction for the backbone dihedral31 was used for the protein. The force field for the spin-labels was taken from the previous work of Sezer et al.32 A weak harmonic restraint with a force constant of 1 (kcal/mol)/Å2 was applied to the backbone of the protein to avoid any large displacement away from the X-ray structure. The Adopted Basis Newton-Raphson minimization was performed for 2000 steps to eliminate initial bad contacts between various atoms. The restrained-ensemble system was simulated in vacuum at a temperature of 300K using Langevin dynamics for 20 ns. The friction was set to 1.0 for all non-hydrogen atoms and an integration time step of 1 fs was used. Bonds involving hydrogen atoms were constrained to their equilibrium values using the SHAKE algorithm.33 For comparison, a conventional (single-copy) MD was also carried out for 20 ns without any restraint. It should be noted that this simulation differs from the restrained-ensemble simulation carried out with a single replica (N = 1) because it does not include any restraint with respect to the ESR/DEER distance histogram data. The nonbonded interactions were smoothly switched off from 10–12 Å using an atom-based cutoff.</p><!><p>It was previously demonstrated that the restrained-ensemble MD simulations converges rigorously toward the result from maximum entropy principle in the limit of strong biasing potential and large number of replicas.24 This is reassuring because it suggests that theoretical constructs built on the restrained-ensemble approach do not corrupt the resulting distribution functions with spurious and unwanted biases, in a maximum entropy sense. As is clear from the above analysis, applications of the maximum entropy principle requires the determination of several coefficients λi, which must be adjusted to match the known experimental data. While the formal task is clear, this optimization process can rapidly become cumbersome and unwieldy, particularly when there is a very large number of data.</p><p>The great advantage of an approach based on restrained-ensemble simulations is to obtain a solution that is formally equivalent to that from maximum entropy principle while avoiding the obligation to determine a large number of unknown Lagrange multipliers. While the computational load is certainly larger in the restrained-ensemble simulations because of the multiple copies, the simulation automatically evolves in a way that is consistent with the maximum entropy principle without any additional intervention by the user. This becomes particularly important in the case of ESR/DEER distance histogram data because there can be a huge body of data to take into account. To make the algorithm more efficient, it was assumed that the different spin-labels were uncorrelated from one another. Such a mean-field approximation is reasonable because the spin-labels in ESR/DEER experiments are typically separated by a sufficiently large distance. To illustrate the method, three histograms are used to restrain the motion of three spin-labels introduced at position 62, 109, and 132 in T4 lysozyme (T4L) using site-directed cysteine mutations. The nitroxide spin-label is shown in Figure 1. T4L is a small soluble protein of 165 residues that has long served as a benchmark system for ESR and SDSL techniques. 5,6</p><p>The T4L system is shown in Figure 2. In the left view is shown the protein with single-copy spin-labels. In the right view is shown a typical snapshot taken from the restrained-ensemble simulations. In the case considered here, the 3 spin-labels are well separated, by distances on the order of 20–40 Å. It is worth noting that, while the samples used in ESR/DEER are fast-frozen from room temperature to liquid nitrogen temperature, it makes sense to consider the system at temperature T because the ensemble of spin-label rotamers is rapidly quenched and the distance histograms reflect the distribution at room temperature. Thus, T4L frozen in this manner faithfully maintain the structural properties it possesses in solution at room temperature. Nevertheless, this issue may require further attention in the future, as there are experimental indications that the freezing process may affect the configurational distribution of the spin-labels.34</p><p>The internal flexibility of the spin-label is of great importance in the interpretation of ESR data in general. Five dihedral angles, χ1 − χ5, connect the nitroxide ring of the spin-label to the protein backbone (Figure 1). The dihedral angle of the disulphide (S-S) bond, χ3, can adopt one of two configuration, 90° and −90°, which are separated by a large energy barrier. A number of X-ray crystallographic and spectroscopic studies aimed at characterizing the rotameric state of the spin label7–13 have suggested that the internal dynamics of the spin-label side chain arise predominantly from the dihedral angles χ4 and χ5. Variable rotameric states were observed for χ1 and χ2 at different sites in T4L, which led to the hypothesis that these dihedrals are hindered by the formation of a hydrogen bond between the sulfur atom of the spin-label and the backbone amide and Cα. 7,13</p><p>Restrained-ensemble MD simulations using 1, 10 and 25 copies for each spin-labels were carried out for 20 ns. The simulation with a single copy corresponds to the standard approach to impose experimental properties. The results from the mean-field restrained-ensemble MD simulations are compared with experimental ESR/DEER data in Figure 3. The restrained simulation converged rapidly toward the experimental histograms. In contrast, the conventional (single copy) unbiased simulations did not accurately reproduce the experimental histograms (this simulation differs from the restrained-ensemble simulation carried out with a single replica with N = 1 because it does not include any restraint with respect to the ESR/DEER histograms).</p><p>To characterize the effect of the number of copies of the spin-labels used in the restrained-ensemble MD simulation, we calculated the rotameric states of the spin-label at the 3 sites in T4L. The distribution of the χ1 and χ2 dihedrals is displayed as scatter plots in Figure 4 for spin-labels at position 62, 109 and 134 in T4L. For all 3 sites in T4L, the number rotameric states visited by the simulations is the smallest with a single replicas and becomes roughly equivalent with 10 or 25 replicas, indicating that the multiple-copy representation has converged for the present system.</p><p>Restraint-ensemble simulations performed with 10 and 25 copies of each spin-label are found to explore more rotameric states along χ1 and χ2 dihedrals when compared to those obtained from conventional single-copy restrained MD simulation. This suggests that the dihedrals are still flexible, even though distance distributions are restrained to experimental distributions in restraint-ensemble simulation. Usually the rotameric states in the spin-label side chains are denoted by p, t, and m, depending on the range of dihedral angles from 0° to 120°, 120° to −120° and −120° to 0°, respectively. spin-label 62 prefer mostly the mm (78%) and mp (18%) rotamers. In contrast, spin-labels 109 and 134 visit multiple rotameric states, with mt (31%), pp (21%), pm (13%), mp (13%), mm (12%) and tm (6%) rotamers for spin-label 109, and mt (42%), mm (25%), pt (14%) and tt (10%) for spin-label 134. The restrained-ensemble MD simulation method enforces the experimental distance histogram while allowing reasonable fluctuations of the dihedral angles of the spin-labels.</p><p>The present treatment implicitly considered the real-space pair-distance ESR/DEER histograms as the "raw input" information from experiments. In reality, such distance histograms are obtained by post-processing the measured time-dependent ESR/DEER signal.35 A legitimate question is whether one should construct the restrained-ensemble simulation methodology on the basis of the actual time-dependent EPR/DEER experimental signal. The observed EPR/DEER signal, which depends on multiple factors, can be expressed as Vij(t) = B(t)[1 − λFij(t)].35 Of primary interest is the function Fij(t), which corresponds to the magnetic interaction between spins i and j. In addition, the constant λ is related to the depth of modulation of the signal and depends on the fraction of doubly labeled protein in the sample. Lastly, there is the background decaying function B(t), which depends on intermolecular spin-spin coupling with distant neighbors. B(t) is sensitive to the dimensionality of the spatial distribution of labeled proteins and is generally unknown. There is a linear relationship between the time-dependent signal, Fij(t), and the distance distribution Hij(r), which can be expressed in discrete form as,</p><p>
(14)Fij(tn)=∑rmM{tn,rm}Hij(rm) where the operator M{tn,rm} includes, among other things, the effect of averaging the spin-spin vector (ri − rj) over all possible orientations. Determining the real-space histograms Hij(r) from the time-dependent signal Fij(t) amounts to a difficult inverse problem, which is typically is solved by using Tikhonov regularization.35 Clearly, converting Vij(t) into Hij(r) requires careful considerations of multiple factors pertaining to the sample preparation and experimental conditions. The validity of a structural model can be ascertained by examining if the time-dependent signal reproduced experimental ESR/DEER data. For example, Boura et al36 generated a large number of possible conformations from unbiased simulations, and then constructed a weighted ensemble by scoring those conformations in post-analysis according to their ability to reproduce the time-dependent ESR/DEER signal. In contrast, the goal of the present multiple-copy restrained-ensemble simulation method is to actively drive the protein conformation toward a 3D structure that is consistent with the data. While a post-analysis scoring method can work if the correct structure is found among the trial set of conformation, to have a practical structural refinement algorithm, it is necessary to have the ability to bias the dynamical evolution of the system toward conformations that are consistent with the input data. To go from ESR/DEER measurements to the 3D protein requires two very different and difficult operations: (1) processing of the time signal, and (2) refining the 3D protein structure from correlated data. While it is conceivable that one could use the signal Vij(t), or the post-processed function Fij(t), directly into some kind of restrained-ensemble MD simulation for the purpose of structural refinement, we see no clear advantage in combining two difficult operations into a single step. In effect, the distance histograms offer a nice intermediate stage allowing one to get a sense of the spatial information provided by the data, and see if the post-processing treatment of the signal has been carried out properly.</p><!><p>A multiple-copy mean-field restrained-ensemble method was implemented to incorporate the information from multiple ESR/DEER distance histograms simultaneously. Illustrative tests on 3 coupled spin-labels in T4L show that the method efficiently imposes the experimental distance distribution rapidly in this system. Future directions of the present methodology will be to fully explore the consequences of incorporating all the ESR/DEER data on the rotamer distribution of the spin-labels in T4L. While the present illustrative simulations were carried out in vacuum, it will certainly be critical to consider the effect of explicit solvent in future studies. It might also be of interest to try to reduce the large sulphur-sulphur dihedral barrier to increase the rotameric interconversion rates of the nitroxide spin-label using enhanced sampling techniques such as accelerated MD (aMD).37 The restrained-ensemble MD simulation methodology described here could be extended to other experimental properties and observables, such as X-ray and neutron scattering profiles, NMR, resonance energy transfer to give a few examples. Future applications will be aimed at determining the population of nitroxide spin-label rotamers using ESR/DEER experiments.</p>