Prediction of Bond Vector Autocorrelation Functions from Larmor

Article pubs.acs.org/JCTC

Prediction of Bond Vector Autocorrelation Functions from Larmor Frequency-Selective Order Parameter Analysis of NMR Relaxation Data Janet S. Anderson,† Griselda Hernández,‡ and David M. LeMaster*,‡ †

Department of Chemistry, Union College, Schenectady, New York 12308, United States Wadsworth Center, New York State Department of Health and Department of Biomedical Sciences, School of Public Health, University at Albany - SUNY, Empire State Plaza, Albany, New York 12201, United States

‡

S Supporting Information *

ABSTRACT: Protein molecular dynamics interpretation of the standard R1, R2, and heteronuclear NOE relaxation measurements has typically been limited to a single S2 order parameter which is often insufficient to characterize the rich content of these NMR experiments. In the absence of exchange linebroadening, an optimized reduced spectral density analysis of these measurements can yield spectral density values at three distinct frequencies. Surprisingly, these three discrete spectral density values have proven to be sufficient for a Larmor frequencyselective order parameter analysis of the 223 methine and methylene H−C bonds of the B3 domain of Protein G (GB3) to accurately back-calculate the entire curve of the corresponding bond vector autocorrelation functions upon which the NMR relaxation behavior depends. The 13C relaxation values calculated from 2 μs of CHARMM36 simulation trajectories yielded the corresponding autocorrelation functions to an average rmsd of 0.44% with only three bond vectors having rmsd errors slightly greater than 1.0%. Similar quality predictions were obtained using the CHARMM22/CMAP, AMBER ff99SB, and AMBER ff99SB-ILDN force fields. Analogous predictions for the backbone 15N relaxation values were 3-fold more accurate. Excluding seven residues for which either experimental data is lacking or previous MD studies have indicated markedly divergent dynamics predictions, the CHARMM36-derived and experimentally derived 15N relaxation values for the remaining 48 amides of GB3 agree to an average of 0.016, 0.010, and 0.020 for the fast limit (Sf2) and Larmor frequency-selective (SH2 and SN2) order parameters, respectively. In contrast, for a substantial fraction of side chain positions, the statistical uncertainties obtained in the relaxation value predictions from each force field were appreciably less than the much larger differences predicted among these force fields, indicating a significant opportunity for experimental NMR relaxation measurements to provide structurally interpretable guidance for further optimizing the prediction of protein dynamics.

1. INTRODUCTION Beyond the specialized applications of time-resolved crystallography,1 no experimental technique is sufficiently robust to simultaneously quantify both the spatial and temporal aspects of protein motion in atomic detail. As a result, conformational dynamics analysis must directly integrate molecular simulation methods as the means by which the relevant experimental data are to be quantitatively interpreted. Conversely, experimental measurements which offer a stringent characterization of protein dynamics behavior can serve to highlight the shortcomings of the computationally practical algorithms currently used to model that motion and thus provide a basis for future optimization of those simulation methods. NMR relaxation measurements offer an unmatched capability for experimentally quantifying key aspects of molecular dynamics with atomic level detail, and in principle, these experimental data are directly predictable from MD trajectories. After correction for possible conformational/chemical exchange linebroadening effects, the standard R1, R2 relaxation © XXXX American Chemical Society

rate and heteronuclear NOE experiments directly monitor the degree of orientational disorder for the individual 1H−15N or 1 H−13C bond vectors. However, quantitative comparison of these experimental results to MD-derived predictions presents significant challenges from both sides of the problem. Inaccuracies in the rotational diffusional behavior of both solute and solvent predicted with the water models in current use are substantially larger than the typical experimental uncertainty for 15N relaxation measurements on small wellbehaved proteins.2 Most commonly, this complication is dealt with by assuming the separability between the global tumbling bond vector autocorrelation functions and the internal motion autocorrelation functions. After superimposing the trajectory frames upon the initial frame and calculating the internal motion autocorrelation functions, a global diffusion tensor is estimated from the experimental relaxation data, and the bond Received: April 11, 2017 Published: May 25, 2017 A

DOI: 10.1021/acs.jctc.7b00387 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

these experiments are directly sensitive to only a small set of frequency components within the overall spectral density functions. Given the fact that the three R1, R2, and heteronuclear NOE values are insufficient to determine the five spectral density values on which they depend, Peng and Wagner10 have described an additional set of NMR relaxation measurements which are formally adequate for a complete determination of the spectral density components. Unfortunately, in practice, these supplementary relaxation experiments are generally less robust than the standard R1, R2, and NOE measurements. To circumvent this limitation, the reduced spectral density analysis11−13 introduces a J(ωH*) value to mimic the averaged effects of the J(ωH-ωX), J(ωH), and J(ωH+ωX) factors. Using microsecond-scale simulations on the G B 3 p r o t e i n w i t h t h e CH A R M M 2 2 / C M A P , 1 4 , 1 5 CHARMM36,16,17 AMBER ff99SB,18 and AMBER ff99SBILDN19 force fields, we have found that the accuracy with which the reduced spectral density analysis back-calculates the MD-derived spectral density values can be significantly improved upon modest alteration of the originally proposed12 values for ωH* of 0.870*ωH and 1.563*ωH in 15N and 13C relaxation studies, respectively. Drawing upon the simplifications offered by the reduced spectral density formalism, this study considers an alternate extension of the Lipari−Szabo paradigm based on a Larmor frequency-selective (Sf2,SH2,SX2) order parameter analysis20

vector orientations from a crystallographic or NMR structural model of the protein are used to predict the individual global tumbling bond vector autocorrelation functions. The products of these individual internal and global motion bond vector autocorrelation functions are then Fourier transformed to yield the spectral density function at the five Larmor frequency combination values of J(0), J(ωX), J(ωH-ωX), J(ωH), and J(ωH+ωX) from which the individual NMR experimental relaxation values can be directly derived, where ωX is the Larmor frequency for either 15N or 13C. Unfortunately, in the case of anisotropic tumbling, inaccuracies in the modeling of internal motion alter the orientation of the bond vectors with respect to the principal axes of diffusion, thus distorting the corresponding global motion bond vector autocorrelation functions as well. Ideally, the experimental NMR analysis attempts to reverse this process by using the R1, R2, and heteronuclear NOE values to estimate the characteristic spectral density values and from those values provide a model for the underlying bond vector autocorrelation functions. In Lipari and Szabo’s original formulation of this analysis,3 a two-parameter (S2,τe) model was proposed in which the order parameter characterizes the degree of orientational disorder that arises from internal bond vector dynamics, while τe is the effective time constant for the formation of that disorder as modulated by the global tumbling process (i.e., 1/τe = 1/τM + 1/τi, where τM and τi characterize global and internal motion, respectively). Unfortunately, when the value of τi approaches τM, the derived S2 values operationally become a function of the global tumbling rate. This complication has led to a number of proposals for approximating the S2 values obtained from MD trajectories as a function of the τM windowing effect.4−6 In a related effect, the derived τi values are often poorly quantified by the experimental data and are typically not reported. As a result, trajectory-based predictions of protein NMR relaxation data are generally compared against only the experimentally determined S2 values. The limitations of such an approach has been elegantly illustrated in a study by Levy and colleagues for which model 15N R1, R2, and NOE values were calculated as a function of S2 and τi for a small isotropically tumbling protein (Table 1).7,8 The range of NMR relaxation

J(ω) = 2/5Sf 2[SH 2SX 2τM /(1 + ω 2τM 2) + (1 − SH 2)τH /(1 + ω 2τH 2) + (1 − SX 2)SH 2τX /(1 + ω 2τX 2)]

in which 1/τX is set to the Larmor frequency of the heteronucleus 15N or 13C and 1/τH is set to the effective 1H Larmor frequency ωH*. The practical utility of this approach to representing the complete spectral density and autocorrelation functions is that the maximal temporal responsiveness to the changes in these functions occurs at the frequencies for which the experimental measurements are most sensitive. An additional conceptual appeal of the Larmor frequency-selective model free analysis is that the linear functional dependence for each of the three order parameters (Sf2,SH2,SX2) allows these parameters to be analytically expressed in terms of the reduced spectral density values which is infeasible for the more widely used (Sf2,Ss2,τs) formalism.9 The molecular dynamics simulation studies demonstrate that the Larmor frequency-selective spectral density function given in eq 1 provides strikingly accurate representations of the bond vector autocorrelation functions derived directly from these simulations for not only the comparatively immobile backbone 1 H−15N bond vectors but also for the considerably more mobile methine and methylene 1H−13C bond vectors. While there is no a priori reason why the autocorrelation functions of these various simulations should precisely conform to this three-parameter equation, individual sets of the experimentally accessible R1, R2, and heteronuclear NOE relaxation values corresponding to each bond vector were sufficient to backpredict the entire original autocorrelation functions to an accuracy of 1% or better in virtually all cases. Appreciably poorer performance was observed for this set of bond vectors when the analogous spectral density equation for the (Sf2,Ss2,τs) formalism was utilized. The practical relevance of this level of performance is enhanced by the fact that methine and methylene groups constitute 89 of the 98 stereochemically distinct protein C−H bond vector types. Unfortunately, to date,

Table 1. Range of Relaxation Values as a Function of the Time Constant for Internal Motion, τia S2 0.80 0.40 a

N R1 (s−1)

15

2.43 to 2.88 1.23 to 2.68

15

N R2 (s−1)

5.4 to 5.9 2.8 to 4.4

15

(1)

N NOE

0.40 to 0.60 −0.60 to 0.30

Here, ωN of 50.7 MHz, τM of 4 ns, τi from 0 to 2 ns.

values compatible with a given S2 value is substantial, particularly as the degree of motional disorder is increased. The original (S2, τe) Lipari−Szabo analysis3 was subsequently extended to a three-parameter formalism consisting of a fast limit order parameter Sf2 and a slower time frame order parameter Ss2 with a corresponding time constant τs.9 Although a wider set of experimental relaxation data could be fitted by this extended formalism, the τs parameter can still be poorly constrained and the physical significance of this internal time constant is often open to question. Ideally, the optimal analysis of NMR relaxation data should provide the most accurate representation of the underlying bond vector autocorrelation functions. On the other hand, B


Article


force fields. The complete autocorrelation functions C(t) were obtained for each bond vector by multiplying Ci(t) by the bond vector-specific global tumbling autocorrelation functions, determined as described below. The values of the spectral density function J(ω) at the requisite frequencies for the NMR relaxation parameter predictions were calculated by summing over the time interval covered by C(t) according to

relatively few NMR relaxation studies have analyzed these bond vector types. Regarding the nine methyl bond vector types (one in alanine, methionine, and threonine; two each for valine, leucine, and isoleucine), the dominant relaxation contribution of methyl group rotation typically occurs at frequencies well above the highest sampled NMR frequency so that the advantages of the Larmor frequency-selective analysis are correspondingly diminished. Given the marked systematic differences in prediction observed among these force fields, particularly for the side chain dynamics, considerable potential can be anticipated for the use of experimental NMR relaxation measurements in guiding the future evolution of force field development.

J(ω) =

(3)

While the digitization artifact correction for the initial time point can be approximated by a 2-fold decrease in the value at t = 0, the alternate approach of estimating the initial increment by a 50% increase in the value at t = 5 ps was employed. This choice reflected the fact that, due to the rapid decay of the highest frequency components, the decrease between C(0 ps) and C(5 ps) can be greater than 30-fold larger than the decrease between C(5 ps) and C(10 ps) for the comparatively restricted bond vectors. For the well-ordered backbone 1H−15N bonds of GB3, the discrepancy between these two approaches to the digitization artifact correction is 0.7% of the total spectral density at ωH*, on the order of the experimental uncertainty for the heteronuclear NOE measurements. In comparing the dynamics modeled by MD simulations with the experimentally observed NMR relaxation values, the theoretical nuclear dipolar coupling factor must be rescaled to reflect the fact that MD simulations standardly utilize rigid models for hydrogen-containing covalent bonds so that the contributions to relaxation arising from high frequency classical/quantum mechanical librations are not adequately predicted.25 While this rescaling of the theoretical coupling factor to optimally match the experimental data is typically expressed in terms of an altered effective 1H−15N bond length, in reality any error in the assumed average 15N chemical shift anisotropy (CSA) value for NMR data at a given field strength is simultaneously corrected for as well.26 Effective bond length optimizations for the 1H−15N data were carried out under the assumption of a uniform 15N CSA value of 168 ppm.24 Under these conditions, an 1H−15N bond length of 1.036 Å was obtained for the CHARMM simulations. 2.3. Rotational Velocity Rescaling. Following translation of each CHARMM36 trajectory frame to the all-heavy-atom center of mass, the superposition rotation between each adjacent pair of frames was determined using the quaternion algorithm described by Kearsley27 as implemented in FORTRAN by Rupp and Parkin (PDBSUP28). The rescaling of the rotational velocity throughout the trajectory was modeled by multiplying the angle of rotation for each step of superposition by a constant factor and then renormalizing the corresponding unit quaternion.24 The fit to the experimental 15 N relaxation data29 was optimized when the rescaling factor was set to 0.644. This change in step size for a Brownian diffusion process alters the resultant global tumbling rate by the square of this factor which compensates for inaccuracies in both the computational modeling of the rotational diffusion behavior and the differences in the effective temperatures for the simulation and experimental data sets. As the motions with frequencies near the global tumbling rate are more weakly sampled than the more rapid dynamics, an exponential extrapolation from 2.5 ns was used to model the tail of the bond vector autocorrelation function.24 2.4. Rotational Diffusion Tensor. The experimental 15N R1, R2, and NOE values for the dynamically restricted backbone

2. METHODS 2.1. MD Simulations. Coordinates from the 1.1 Å resolution X-ray structure (PDB code 2IGD21) were modified with CHIMERA22 to form the Δ1−5, T6M, T7Q variant that has been widely used for NMR studies. Hydrogen atoms were added to the protein heavy atoms and crystallographically defined water molecules with all carboxyl and aliphatic amine groups being charged. The protein coordinates were placed in a rectangular box with no protein atoms within 10 Å of the boundary. The simulation box was filled with TIP3P waters (modified version for CHARMM simulations). The starting configuration for each trajectory was then generated by an independent placement of 11 sodium ions and 9 chloride ions added to establish electroneutrality and an ionic strength near 150 mM. The simulations of GB3 were carried out using CHARMM22/CMAP and CHARMM36 force fields (2 μs total for each) in the NAMD2 program23 and using AMBER ff99SB and AMBER ff99SB-ILDN force fields (1.5 μs total for each) in the AMBER16 program. Particle mesh Ewald summation was applied for long-range electrostatics using a 1 Å grid spacing. Short range interactions utilized a switching distance of 10 Å and a cutoff of 12 Å. SHAKE constraints were applied for all bonds to hydrogen, and a step size of 1 fs was used. The protocol for equilibration at 298 K was carried out as previously described.24 The production phase was carried out under NVE conditions with the first 10 ns being treated as equilibration. Production runs of 500 ns were carried out with each force field excepting the AMBER ff99SB-ILDN force field for which four modestly shorter runs were combined for a total of 1.5 μs of simulation time. Trajectory frames were stored at 5 ps intervals. For the CHARMM simulations, the temperature drift over the trajectory was minimal with the averages over consecutive 10 ns intervals remaining within 0.1 degrees of the overall average temperature. Larger monotonic drifts in temperature of a few degrees were observed in the AMBER simulations. 2.2. Bond Vector Autocorrelation and Spectral Density Functions. Following superposition of the trajectory frames upon the initial frame, the H−N and H−C bond vector autocorrelation functions were calculated according to Ci(t ) = 3[(e(τ )(e(τ + t )]2 − 1 /2

∑ 2C(t )cos(ωt )Δt

(2)

averaged over each trajectory where e(t) is the unit vector defining the bond orientation. The autocorrelation functions were calculated out to 20 ns (∼6-fold larger than the global tumbling time). A measure of the quality of statistical sampling was indicated by the fact that no deviations from monotonicity in the Ci(t) values greater than 0.002 were observed for any of the amide, methine, or methylene bond vectors with any of the C


Article


τM). The complete autocorrelation functions were then approximated by the product of the MD-derived internal autocorrelation functions and the global autocorrelation functions derived from the experimentally determined diffusional ellipsoid. Discrete Fourier analysis yielded the corresponding spectral density values J(ω) which were then used to calculate the R1, R2, and NOE values arising from the 1H−15N dipolar and 15N CSA interactions36

amides of GB3 were analyzed with the ModelFree 4.1 software30 under the assumption of a prolate rotational diffusion ellipsoid using the 1.1 Å resolution crystal structure (PDB code 2IGD21) to specify the angle between each bond vector and the principal axis of the diffusion tensor. To predict the contribution to orientational disorder arising from global molecular tumbling, the angle of each bond vector with respect to the principal axis was determined for every sampled trajectory frame, and the relative instantaneous rotational diffusion rate was predicted from the initial slope of the Woessner equation for axially symmetric diffusion31,32 C(t ) =

⎛ d2 ⎞ R1 = ⎜ ⎟[J(ωH − ωX ) + 3J(ωX ) + 6J(ωH + ωX )] ⎝4⎠

1 ⎡⎣(3cos2θ − 1)2 exp[− 6Dx t ] 4 2

+

2

(5)

⎛ d2 ⎞ R 2 = ⎜ ⎟[4J(0) + J(ωH − ωX ) + 3J(ωX ) + 6J(ωH ) ⎝8⎠

+ 12cos θ sin θ exp[− (5Dx + Dz )t ] + 3sin 4 θ exp[− (2Dx + 4Dz )t ]⎤⎦

(Δσ ·ωX )2 [J(ωX )] 3

(4)

These relative instantaneous diffusion rate estimates were averaged across the trajectory and then used to identify the corresponding effective angle with respect to the principal tensor axis. These angles were then applied to the Woessner equation to model the bond vector-specific global tumbling autocorrelation functions. The experimentally observed isotropically averaged global tumbling time was divided by the relative averaged instantaneous rotational diffusion rates to obtain the bond vector-specific τM values used in the order parameter analysis.

+ 6J(ωH + ωX )] +

(Δσ ·ωX )2 [4J(0) + 3J(ωX )] 18 (6)

NOE = 1 +

d 2 γH · [6J(ωH + ωX ) − J(ωH − ωX )] 4R1 γX

(7)

−3

where d = (μo/4π)ℏγHγX⟨rHX ⟩, γ is the magnetogyric ratio, and Δσ equals the anisotropy of the chemical shift tensor (σ∥ − σ⊥). The CHARMM36-derived relaxation values for 37 residues drawn from the more highly ordered segments of the protein were found to yield quite satisfactory fits to the experimental R1, R2, and NOE values (Table 2).The analogous predictions

3. RESULTS AND DISCUSSION 3.1. NMR Relaxation Predictions for GB3 and the Separability of Internal and Global Dynamics. Most backbone amides of GB3 exhibit highly restricted internal dynamics within the time frame of global tumbling, and the corresponding experimental 15N S2 order parameters are reasonably accurately reproduced by a number of the current force fields.4,33,34 Three internal segments of the chain exhibit more extensive motion on the ∼ns time frame for which various force fields have yielded markedly differing predictions of the dynamics which in many cases do not closely match the experimental S2 values. These regions of discrepancy have been ascribed, in part, to the rearrangement of specific local interactions that occur on a time frame slower than molecular tumbling.4,34 Under the assumption of separability between the internal and global motions, the experimental 15N R1, R2, and NOE values of Bax and colleagues29 were analyzed with the ModelFree 4.1 software30 using the reported 1.1 Å resolution crystal structure (PDB code 2IGD21) to yield an isotropically averaged rotational correlation time of 3.115 ns for a prolate diffusional ellipsoid with an asymmetry ratio (Dz/Dx) of 1.35 from which a corresponding global tumbling autocorrelation function for each bond vector can be directly calculated.31 TIP3P-solvated simulations of GB3 were carried out using CHARMM22/CMAP and CHARMM36 force fields (2 μs total for each) in the NAMD2 program23 and using AMBER ff99SB and AMBER ff99SB-ILDN force fields (1.5 μs total for each) in the AMBER16 program. This length of simulation was sufficient to provide the greater than 100-fold sampling needed to properly evaluate the relatively high frequency motions which occur within the time frame of molecular tumbling.35 After superimposing the trajectory frames upon the initial frame, the various bond vector internal autocorrelation functions were calculated out to 20 ns (>6-fold larger than

Table 2. Force Field Dependence of GB3 15N Relaxation Predictions Averaged over 37 Motionally Restricted Residues

a

force field

Δ15N R1 (%)

Δ15N R2 (%)

Δ15N NOE

CHARMM36 CHARMM22/CMAP ff99sb-ILDN ff99sb CHARMM36-rvra

1.47 1.40 2.25 2.33 1.70

2.00 1.95 2.35 2.78 3.02

0.014 0.009 0.016 0.021 0.024

Rotational velocity rescaled trajectory.

obtained using the CHARMM22/CMAP yielded comparable rmsd values for the fits to the experimental relaxation values, while those derived from the AMBER force fields were somewhat larger. This quality of fits between the predicted and observed relaxation values strongly suggests the absence of any appreciable contribution to the R2 values arising from exchange linebroadening dynamics as has been more directly demonstrated using magnetic field-dependent relaxation measurements.26 Previous model calculations have offered support for the reasonableness of assuming an approximate separability between the internal and global molecular motion for the backbone resonances of well-ordered globular proteins.32 The present simulations can be used to provide a more stringent test for the upper bound of errors arising from the assumption of separability. Recently, we have introduced a rotational velocity rescaling protocol to compensate for the inaccurate diffusional properties of the standard MD water models and the resultant alteration in global molecular tumbling dynamics.24 D


Article

Journal of Chemical Theory and Computation The stepwise global reorientation between adjacent sampled frames of the trajectory is characterized by the quaternion specifying the axis of rotation and angle of rotation around that axis. The angles of rotation throughout the trajectory are then uniformly rescaled in magnitude and used to “rewind” the trajectory so as to match the experimental global rotational correlation time while at the same time preserving any correlation that might have occurred between alterations in molecular conformation and the instantaneous rotational velocity of global tumbling. The CHARMM36 trajectories were rotationally rescaled so as to provide an optimal match to the relaxation values for the same set of 37 highly ordered residues. The rmsd fits to the experimental R1, R2, and NOE values (Table 2) were somewhat larger than for the values obtained using a best fit diffusional ellipsoid and the internal autocorrelation functions obtained from the CHARMM36 simulations, although they were comparable to the quality of fits obtained with the two AMBER-based force fields. As could be anticipated from the more limited statistical sampling of the lower frequencies characteristic of global tumbling, the predicted R2 values exhibited comparatively larger errors in the CHARMM36-rvr analysis (Table 2). A similar conclusion is obtained when the predicted NMR relaxation values obtained from CHARMM36 trajectories were compared between the rotational velocity rescaled analysis and those obtained by superposition and multiplication of the internal and global tumbling autocorrelation functions. The rmsd fits of 0.69%, 1.38%, and 0.009 for the R1, R2, and NOE values, respectively, are less than the deviations from the experimental values and exhibit an increased uncertainty for predicted R2 values. The analogous rmsd values between the NMR relaxation predictions from the superimposed and the “rewound” trajectories for the 223 methine and methylene bond vectors in the CHARMM36 simulations are 1.8%, 2.5%, and 0.007, respectively. Since the occurrence of correlation between the internal motion and the instantaneous rotational diffusion should be preserved in the velocity rescaling analysis, these results indicate that the assumption of separability for the internal and global motion of GB3 yields errors which appear to be less than other sources of force field-dependent error. The relaxation values for Gly 14 derived from the CHARMM-based simulations (Figure 1) most notably differ from the both experimental results and the AMBER-based predictions (Figure S1) in predicting a strikingly elevated degree of conformational dynamics at this position. In contrast, the CHARMM36 trajectories provide more accurate relaxation predictions for Gly 41, while both AMBER-based force fields yield exaggerated estimates of internal dynamics for this residue. The erroneously predicted dynamics at Gly 9, Ala 20, and Asn 22, that have previously been structurally analyzed for the AMBER ff99SB, ff99SB-ILDN, and OPLS-AA force fields,4,34 are also not well predicted in the CHARMM trajectories and were suggested to reflect inadequacies in representing hydrogen bonding interactions with these nonpolarizable force fields. In contrast to the backbone dynamics for which both experiment and simulation exhibit substantial mobility on the ps−ns time frame for only small segments of the protein, the great majority of side chain positions predict extensive dynamics. Encouragement that the present simulations might provide a useful level of conformational sampling for the side chain dynamics can be drawn from 1 ms MD simulation of bovine pancreatic trypsin inhibitor (BPTI) by Shaw and

Figure 1. 15N R1, R2, and NOE values for GB3 at 600 MHz 1H were predicted from 2 μs of simulations with CHARMM22/CMAP (■) and CHARMM36 (⧫) force fields in the NAMD2 program. The experimental relaxation values29 are denoted as circles. The ModelFree 4.1 program30 was used to estimate a prolate diffusion tensor with an isotropically averaged τM of 3.115 ns and an asymmetry of 1.35 from the residues that were adequately fitted by an S2 only dynamical model.

colleagues.37 They identified a set of five distinct backbone conformational basins which slowly interchanged with time constants in the range of μs or longer. Within these individual backbone conformational basins, the average dynamical content for the side chains dropped to near zero for autocorrelation function times above ∼20 ns. Using the rotational diffusion ellipsoid derived from the 15N relaxation measurements and an analogous rescaling of the effective H−C bond length, the 13C relaxation values for all such bond vectors were predicted from each of the four force fields under study. The enhanced dynamics of side chains is illustrated for each methine or methylene bond vector for the Cβ positions in the AMBER ff99SB and ff99SB-ILDN trajectories (Figure 2). This comparison is particularly informative given the fact that the force field reparameterization in the ff99SB-ILDN analysis19 was limited to adjusting the side chain torsion angle potentials for the isoleucine, leucine, aspartate, and asparagine residues. For the majority of residues, the independent simulation studies using the two AMBER force fields yielded quite similar side chain NMR relaxation predictions. On the other hand, a subset of residues exhibit quite substantial differences between the two force field predictions. The majority of those residues yielding substantial differential dynamics predictions belong to the four reparameterized residue types. It should be noted that for a number of the methylenebearing residues, the two geminal H−C vectors exhibit appreciably different relaxation values. Such nonequivalence can reflect the presence of asymmetry in the local conformaE


Article


tumbling rate. The reverse holds for the bond vectors whose reference orientation lies in the plane of the two minor axes of the tensor. To account for this effect, the effective angle for each bond vector was determined by averaging the instantaneous relative diffusion rates predicted from the Woessner diffusion equations31,32 applied to bond vector angles found in each frame of the superimposed trajectory. Insight into the operational significance of these effective bond vector orientations can be gained by considering the methine and methylene bond vectors of GB3 which lie closest to the principal axis (25% of total) or nearly perpendicular to that axis (25% of total) in the high resolution crystal structure. The differences in the bond vector-specific rotational diffusion times from the isotropically averaged τM value were compared for the trajectory-averaged instantaneous diffusion times vs those derived from the crystallographically determined orientations (Figure 3). As the degree of orientational disorder increases,

Figure 2. 13C R1, R2, and NOE values (600 MHz 1H) for the 1 β 13 β H − C methine and methylene bonds of GB3 were predicted from 1.5 μs of simulations with the AMBER ff99SB (■) and ff99SB-ILDN (⧫) force fields in the AMBER16 program. The rotational diffusion tensor determined from the 15N relaxation measurements and crystal structure21 was applied. As the 13C CSA values for these positions are comparatively small, only dipolar interactions were included using an effective 1H−13C bond length of 1.103 Å.

tional dynamics as monitored by these two approximately orthogonal bond vectors, as has previously been noted in relaxation studies of the geminal methyls of leucine residues.38,39 On the other hand, these nonequivalences may also reflect differences in the orientation of these bond vectors with respect to the rotational diffusion tensor. Applying the assumption of separability between global and internal motion for an anisotropically tumbling protein operationally requires specifying an effective orientation for each bond vector with respect to the principal axes of the rotational diffusion tensor so as to derive a global tumbling autocorrelation function for each bond vector. Current NMR relaxation analysis software typically derives the diffusion tensor utilizing a high resolution structure to define the orientation of each bond vector relative to the principal axes of the tensor. For a protein such as GB3 with a rotational asymmetry of 1.35, the effective tumbling rates for different bond vectors can differ by up to 17.5% with the corresponding R1 and R2 relaxation values varying near proportionally for the motionally restricted bond vectors. For a position exhibiting significant internal mobility, the use of crystallographically determined bond orientations can be problematic. Given the conceptual utility of assuming axial symmetry during the analysis of internal molecular motion,3 determining the simple time-averaged bond orientation within the frames of the superimposed trajectory might initially seem appropriate for defining the effective orientation with respect to the diffusion tensor. However, for a bond that is oriented along the principal axis of a prolate diffusion tensor, all reorientations of that vector will give rise to a more rapid

Figure 3. Averaging of 1H−13C bond vector internal orientation during global rotational diffusion. The methine and methylene bond vectors of GB3 most closely oriented parallel (25%) or perpendicular (25%) to the principal axis of the prolate diffusion tensor in the crystallographic structure were analyzed. The effective angles with respect to the principal axis were determined for each vector by averaging the relative instantaneous rotational diffusion rates predicted for each trajectory frame. The difference in the bond vector-specific τM values from the isotropically averaged τM were normalized to the analogous difference obtained from the crystallographically determined orientations ([τM(MD) − τM(iso)]/[τM(X-ray) − τM(iso)]) and compared to the aggregate order parameters for that bond vector.

the trajectory-averaged values generally converge toward the isotropic value. These results emphasize the fact that inaccuracies in the force field modeling of protein internal motion will enter into the prediction of experimental NMR relaxation values not only in the analysis of the internal dynamics but also in the representation of global tumbling for anisotropic systems. 3.2. Reduced Spectral Density Analysis. The three high frequencies for which the spectral density function is sampled during 15N relaxation (ωH−ωX, ωH, and ωH+ωX) differ by a factor of only ±10%. The reduced spectral density approximation sets these high frequency spectral values to a common F


Article

Journal of Chemical Theory and Computation effective frequency value which then allows the remaining spectral density values to be derived from experimental relaxation measurements11−13 J(0) = [R 2 − 0.5R1 − 0.6 ·RNOE]/[d 2/2 + 2(Δσ ·ωN )2 /9] (8)

J(ωN ) = [R1 − 1.4·RNOE]/[3d 2/4 + (Δσ ·ωN )2 /3]

J(ωH *) = RNOE /[(5d 2/4)]

(9) (10)

where RNOE = (γN/γH)R1(NOE − 1) Under the assumption that the spectral density function takes the form given by J(ω) = λ1/ω2 + λ2, Kay and colleagues12 demonstrated that the NOE factor [6J(ωH+ωN) − J(ωH−ωN)] should be best approximated with a ωH* value of 0.870ωH. While differing values for an effective average ωH* are obtained from the equations for R1 and R2 (eqs 5 and 6), these relaxation values are appreciably less dependent on the high frequency spectral densities. A similar analysis of 13C relaxation contributions yields an NOE-optimized effective ωH* value of 1.563ωH. While a number of studies have considered how robustly the assumption of a single effective 1H frequency enables protein 15N relaxation analysis,40,41 the analogous effects in 13C relaxation studies have been less extensively examined.42,43 An assessment of the practical performance of the reduced spectral density analysis is readily performed by calculating the NMR relaxation values (eq 5−7) from the spectral density functions determined directly from the simulation trajectories, using these relaxation values to determine the reduced spectral density values (eq 8−10), and then comparing the derived J(ωH*) to the corresponding spectral density value obtained from the initial trajectory analysis. When the value of 0.870ωH was assumed for ωH*, the derived J(ωH*)RSD values were systematically 2% higher than those obtained from the original GB3 trajectories. Adjusting ωH* to 0.879ωH yielded a unit slope with an rmsd of 0.9% for the ratios of corresponding J(0.879ωH) values for the 55 backbone amides in the CHARMM36 (Figure 4). Highly similarly qualities of fit were obtained for the J(0.879ωH) values using each of the other three force fields examined. While, in principle, the J(0) and J(ωN) values derived from the reduced spectral density analysis could also deviate from the corresponding spectral density values obtained from the initial trajectory analysis, their agreement was more precise than that obtained for the J(0.879ωH) analysis. Application of the reduced spectral density analysis to 13C relaxation data can be anticipated to be more problematic. In this case, both ωH−ωX and ωH+ωX differ from ωH by 25%, a factor which enters into the effective frequency estimate approximately as the square. When the J(ωH*) analysis was carried out for the 223 methine and methylene bond vectors of the CHARMM36 trajectories assuming an optimal ωH* equal to 1.563ωH,12 the derived J(ωH*)RSD values again systematically deviated from the original simulation values. A unit slope was obtained when ωH* was set to 1.541ωH with an rmsd of 3.9% for the ratios of corresponding J(1.541ωH) values (Figure 5A). While obviously an appreciably greater variability than that obtained for the 15N analysis, this dispersion is over 70-fold less than [(ωH+ωC)/[(ωH−ωC)]2 which approximates the range of values being averaged within each bond vector spectral density function. The corresponding analysis of the CHARMM22/ CMAP trajectories also yielded a unit slope for 1.541ωH. For

Figure 4. Correlation between the spectral densities at a frequency equal to 0.879 times the 1H Larmor frequency for the backbone amide bond vectors of GB3 determined directly from the CHARMM36 trajectories and determined from applying the reduced spectral density analysis to the 15N R1, R2, and NOE values predicted from the CHARMM36 trajectories.

the AMBER trajectories the unit slope was observed at 1.539ωH, although the effect of such a modest shift on the overall results appears to be negligible. While the 13C reduced spectral density analysis of the CHARMM36 trajectories yielded J(0)RSD values that were closely matched to those derived from the original trajectories, the J(ωC)RSD values were systematically 7% too high. This deviation reflects the fact that the ωH* value of 1.541ωH being used is derived from the 6J(ωH+ωX)−J(ωH+ωX) dependence of the NOE formula (eq 7). The corresponding predicted ωH* values in the 13C reduced spectral analysis are 1.058ωH and 1.116ωH when optimized for J(0) and J(ωC), respectively.12 In the Larmor frequency-selective order parameter analysis described below, a 7% elevation in the initially predicted J(ωC)RSD values yields only a modest systematic deviation in the estimated order parameters. To exploit the utility of this initial prediction, the first round of order parameter predictions was then used to re-estimate the spectral density functions with the derived values at 1.058ωH and 1.116ωH being substituted for ωH* in the iterated calculations of J(0) and J(ωC), respectively. This iterated reduced spectral density analysis of 13C relaxation decreased the rmsd on the J(0)RSD/J(0)MD values more than 2-fold to a level of 0.7% over the 223 methine and methylene bond vectors. More importantly, the average elevation in the iteratively predicted J(ωC)RSD values was reduced to 1.4% with an rmsd of 3.1% (Figure 5B). The plot of J(ωC)RSD vs J(ωC)MD indicates two major populations having distinct slopes. For the lower-sloped population, the J(ωC)RSD and J(ωC)MD values are quite precisely equal. For the indicated set of bond vectors, all of the predicted NOE values are in the rigid tumbling limit (≤1.19) with the variation of J(ωC) values principally reflecting the bond vector orientations with respect to the rotational diffusion tensor. For most bond vectors exhibiting NOE values outside of the rigid tumbling limit, an additional deviation arises which is approximately proportional G


Article


3.3. Larmor Frequency-Selective Order Parameter Analysis of Spectral Density Functions. Presentation of the Larmor frequency-selective spectral density function (eq 1) can be appreciably simplified by denoting the (1 + ω2τ2) factors as Dab, where a and b denote the frequency and time components, respectively. So that J(0) = 2/5Sf 2[SH 2SX 2τM + (1 − SH 2)τH + (1 − SX 2)SH 2τX ] (11)

J(ωH ) = 2/5Sf 2[SH 2SX 2τM /DHM + (1 − SH 2)τH /DHH + (1 − SX 2)SH 2τX /DHX ]

(12)

J(ωX ) = 2/5Sf 2[SH 2SX 2τM /DXM + (1 − SH 2)τH /DXH + (1 − SX 2)SH 2τX /DXX ]

(13) 2

Adding the further denotation of J′(ω) = J(ω)/(2/5Sf ), these equations can be combined to eliminate the first and third terms. The remaining (1 - SH2)τH term yields [J ′(0) − DHM J ′(ωH )]/(1 − DHM /DHX ) − [J ′(0) − DXM J ′(ωX )]/(1 − DXM /DXX ) = [(1 − DHM /DHH )/(1 − DHM /DHX ) − (1 − DXM /DXH )/(1 − DXM /DXX )](1 − SH 2)τH (14)

and SX 2 = [J ′(0) − (1 − SH 2)τH − SH 2τX ]/[SH 2(τM − τX )] (15) 2

As Sf functions as a scale factor on the spectral densities, the SH2 and SX2 values are determined as a function of the Sf2 value. While fitting on J(0) in eq 11 is largely sufficient for optimizing the Sf2 values, somewhat better performance in prediction of the trajectory-based autocorrelation functions was obtained when the rmsd errors for fitting J(ωH) and J(ωX) values from eq 12 and 13 were included. While the frequency values at which the spectral density function is to be sampled are determined by the reduced spectral density analysis formalism and the global tumbling time τM is experimentally specified, in principle, the τH and τX parameters can be set arbitrarily. The Larmor frequencyselective order parameter analysis was proposed on the premise that the optimal utility in representing the underlying spectral density/autocorrelation functions can be achieved by setting these two constants to 1/ωH* and 1/ωX, respectively, so as to maximize the parameter responsiveness for the frequencies at which the NMR experiments are most sensitive. Under that assumption, with a ωH* value of 0.879ωH for the 15N analysis, DHX and DXH equal 76.18 and 1.0133, respectively, and ωH*/ ωN is 8.67. Using 1.541ωH for the 13C analysis yields DHX and DXH values of 38.57 and 1.0266, respectively, and ωH*/ωC is 6.13. In each case, Daa is simply equal to 2. To examine the robustness of this premise, the order parameters derived for the 223 methine and methylene bond vectors were used to predict the original trajectory autocorrelation functions as determined assuming separability between the internal and global tumbling dynamics. As summarized in Table 3, the bond vector autocorrelation functions from the CHARMM36 simulations were fitted over the interval from 1.00 to 0.01, using the derived Sf2, SH2, and

Figure 5. Correlation between the calculated and back-predicted spectral densities for the methine and methylene 1H−13C bond vectors of GB3. In panel A, the spectral densities were calculated at a frequency equal to 1.541 times the 1H Larmor frequency determined directly from the CHARMM36 trajectories and determined from applying the reduced spectral density analysis to the 13C R1, R2, and NOE values predicted from the CHARMM36 trajectories. In panel B, the spectral densities were calculated at the 13C Larmor frequency determined directly from the CHARMM36 trajectories and determined from applying the reduced spectral density analysis to the 13C R1, R2, and NOE values predicted from the CHARMM36 trajectories, following an iteration for the J(0) and J(ωC) values using the effective 1H frequencies optimized for the R1 and R2 predictions. The brackets indicate 1Hα−13Cα bond vectors with NOE values near the rigid limit (1%

−Lys Hε

0.44 0.51 0.42 0.45 0.51

0.48 0.55 0.50 0.52 0.57

1.42 1.23 2.90 1.56 1.52

3 9 7 11 14

1 9 3 9 12

Here, 223 methine and methylene vectors with fitting of C(t) from 1.0 to 0.01. a

SC2 order parameters, to yield a mean rmsd value of 0.44% when τH is set to 1/ωH* and τC is 1/ωC. Only three H−C vectors yielded rmsd fits to the original autocorrelation functions that were above 1%, the largest two values arising from the Hε−Cε bonds of the highly dynamic Lys 19 side chain. The Hβ2−Cβ bond of Asn 35 yielded the next worst fit with an rmsd of 1.04% (Figure 6A). In contrast, when the analogous fit to the reduced spectral density values J(0), J(ωC), and J(ωH*) were carried out using the more widely studied (Sf2,Ss2,τs) formalism9 (grid search intervals of 0.001 for Sf2 and Ss2, 1 ps for τs), 46 of the 209 methine and methylene bond vectors not involving lysine Hε−Cε bonds gave rmsd fits to the original autocorrelation functions that were above 1%. To illustrate the effects of sampling statistics on the Larmor frequency order parameter modeling of the simulation trajectories, the four 500 ns CHARMM36 trajectories used in this study were analyzed individually for the dynamics of the comparatively ill-fitted Hβ2-Cβ bond of Asn 35. The dispersion among these back-predicted autocorrelation functions (Figure 6B) were comparable to the discrepancy found between the modeled and the trajectory-derived autocorrelation functions for the full 2 μs simulation period (Figure 6A), consistent with the magnitude of the fitting error in this order parameter analysis being similar to the reproducibility between separate trajectory runs. To examine whether setting the two time constants to 1/ωH* and 1/ωC yields an optimal performance in back-predicting the original autocorrelation functions, analyses were conducted for altered time constant values. Decreasing the τH value by 20% provided improved fits for the highly dynamic lysine Hε−Cε bonds at the cost of an increase both in the mean rmsd value and in the number of substantial outliers among the full set of methine and methylene bond vectors. The other analogous shifts of the two time parameters also tend to degrade the performance of the autocorrelation function predictions (Table 3). While the quality of the fits are marginally worse when the corresponding analysis was applied to the ff99SB-ILDN trajectory results (Table S1), a similar conclusion can be drawn that the setting of τH to 1/ωH* and τC to 1/ωC provides a most robust choice for optimizing the utility of the Larmor frequency-selective spectral density representation. Note that for applying these same order parameters to the analysis of the superimposed trajectories, the effective internal time constant 1/τHi = 1/τH − 1/τM must be used and similarly for τXi.3 The (Sf2,SH2,SC2) order parameter analysis for mobility at the Cβ positions in the CHARMM simulations indicates a broad complexity of conformational dynamics (Figure 7). In the fastest time frame, as represented by the Sf2 order parameter (< ∼35 ps), the conformational dynamics at these positions is rather uniformly restricted in angular dispersion. A subset of

Figure 6. Back-predicted autocorrelation functions derived from the 1 β2 13 β H − C bond vector for Asn 35. Excepting the highly mobile 1 ε 13 ε H − C bonds of Lys 19, this is the poorest fit obtained among the 223 methine and methylene bond vectors of GB3. The complete autocorrelation function predicted from the CHARMM36 trajectories and the experimentally determined rotational diffusion tensor, sampled at 5 ps intervals, is shown in red (panel A). The global tumbling autocorrelation function for that bond vector is shown in black. In blue is indicated the back-predicted autocorrelation function obtained by applying the reduced spectral density method to the initially calculated NMR relaxation values, followed by Larmor frequency-selective order parameter analysis. In panel B are illustrated the autocorrelation functions for this bond vector derived from four sets of the Larmor frequency-selective order parameters predicted from the separate 500 ns CHARMM36 trajectories used in the overall 2 μs trajectory analysis.

side chains begins to exhibit substantial angular sampling in the time frame of τH (172 ps for 600 MHz 1H), while the large majority undergo extensive conformational sampling in the time frame of τC (1.055 ns). Although similar on this qualitative level, the CHARMM22/CMAP and CHARMM36 simulations yield substantially different order parameter predictions for many of the residues. Furthermore, in a significant proportion of those instances, the uncertainty estimated from the dispersion among predictions derived from the individual 500 ns trajectories are appreciably smaller than the difference between the two force field predictions, indicating the potential I


Article


the side chain torsion potentials of those four residue types (Figure 8). In this case, the two force fields provided quite

Figure 7. Larmor frequency-selective order parameters (Sf2,SH2,SC2) for the 1Hβ−13Cβ methine and methylene bond vectors of GB3 predicted from the CHARMM22/CMAP (■) and CHARMM36 (⧫) trajectories. Here, S2 is defined as the product of the three order parameters. Uncertainty estimates were derived from the dispersion among the analogous calculations on the four individual 500 ns trajectories used in the overall analysis.

Figure 8. Larmor frequency-selective order parameters (Sf2,SH2,SC2) for the 1Hβ−13Cβ methine and methylene bond vectors of GB3 predicted from the AMBER ff99SB (■) and AMBER ff99SB-ILDN (⧫) trajectories.

utility of this approach for distinguishing differences in the dynamic behavior of closely related force fields. A note should be made of the dramatically larger uncertainty intervals for the order parameter predictions at Lys 31 and Asn 35 for the CHARMM22/CMAP trajectories which arose from anomalous behavior in one of the four 500 ns trajectory runs. Shortly after the beginning of this run, the side chain of Asn 35 underwent a rotamer transition to the typically rare g+ χ1 torsion angle such that its amide group now formed a hydrogen bond to the carbonyl oxygen of Val 39. The Asn 35 side chain remained in this non-native conformation throughout the rest of the trajectory. When the χ1 torsion angle of Lys 31 was in the trans rotamer state, the aliphatic portion of that side chain contacted the Cβ of the Asn 35 side chain. This interaction appears to account for the alteration in the g−/t rotamer population ratio for the Lys 31 side chain in this anomalous trajectory. The ability to distinguish differences in dynamics behavior among closely related force fields is even more clearly illustrated in the comparison of the predicted order parameters at the Cβ positions in the AMBER ff99SB and ff99SB-ILDN force fields which differ from each other only with regards to

similar order parameter predictions for the majority of residues. The residues exhibiting substantially differing order parameter predictions between the two force fields were primarily either from among modified isoleucine, leucine, aspartic acid, and asparagine side chains or immediately adjacent to such sites. In particular, only Asp 22, Asp 36, and Asp 40 predicted ΔSf2 values greater than 0.05. Leu 12, Lys 13, Asp 22, Thr 25, and Asp 36 predicted ΔSH2 values greater than 0.10, while Thr 11, Leu 12, Lys 13, Thr 17, Val 21, Asp 22, Thr 25, Asn 35, Asp 36, Asn 37, Val 39, and Asp 40 predicted ΔSC2 values greater than 0.10. 3.4. Corrections to the Effective 15N Time Constant Used in Analysis of Rapidly Tumbling Proteins. The initial back-prediction of the CHARMM36 autocorrelation functions from the derived (Sf2,SH2,SN2) order parameters yielded errors near 1% for the highly mobile Gly 14 and Gly 41 residues. This level of discrepancy seemed surprising, given the robustness of the reduced spectral density analysis of the 15N relaxation values. It quickly became apparent that the principal limitation was related to the τN factor. At a magnetic field strength of 14.1T (600 MHz 1H), 1/ωN is 2.62 ns, which is only 15% less than the experimentally measured global J


Article


3.5. Crystallographic- vs Simulation-Based Bond Vector Orientations in Analysis of Experimental NMR Relaxation Data. (Sf2,SH2,SN2) order parameters calculated from the NMR relaxation values derived from the CHARMM36 trajectories were compared against those calculated from experimental NMR data.29 In contrast to the preceding model calculations of this study in which the rotational diffusion tensor and the orientation of each bond vector with respect to that tensor are treated as given, the reliability with which experimental data can be analyzed will depend upon the accuracy of that diffusion tensor and the associated internal vector orientations as well as the robustness of the separability assumption. The order parameters predicted from the reported GB3 15N relaxation values were analyzed using rotational diffusion rate values for each bond vector assuming either the X-ray derived bond orientations or those obtained by averaging the instantaneous diffusion rates over the CHARMM36 trajectories (Figure 10). Each approach to estimating the effective

correlation time for GB3. As is well recognized, robust determination of a time constant for internal motion becomes increasingly problematic when that time constant approaches the time frame of global tumbling. To identify an optimal tradeoff in the Larmor frequency-selective order parameter analysis, the time constant perturbation approach described above for the 13C NMR relaxation analysis was applied to the 15N data. As the τN value used in the order parameter calculations was decreased, both the mean and standard deviation of the fits to the 55 backbone 1H−15N bond vector autocorrelation functions decreased until they reached a minimum at τN = 0.59/ωN where the mean of the fits to the autocorrelation functions was 0.13% and the standard deviation was 0.15%, roughly 3-fold more precise than the corresponding 1H−13C bond vector analysis. The optimal fits to the highly mobile Gly 14 and Gly 41 were reached at slightly lower and slightly higher τN values, respectively. This optimum at τN = 0.59/ωN corresponds quite closely to τM/2. Under these conditions, the order parameters derived for Gly 14 provide the poorest fit to the original autocorrelation function (rmsd of 0.57%) for any of the backbone amides in the CHARMM36 simulations (Figure 9).

Figure 9. Back-predicted autocorrelation functions derived from the most poorly fitted 1H−15N bond vector within the CHARMM36 analysis (Gly 14). The global tumbling autocorrelation function for that bond vector derived from the NMR relaxation measurements29 and crystal structure21 is shown in black. Multiplication by the trajectory-derived internal motion autocorrelation function yields complete autocorrelation function (red), as sampled at 5 ps intervals. In blue is indicated the back-predicted autocorrelation function obtained by applying the reduced spectral density method to the initially calculated NMR relaxation values, followed by Larmor frequency-selective order parameter analysis.

Figure 10. Larmor frequency-selective order parameters (Sf2,SH2,SN2) for the 1H−15N bond vectors of GB3 predicted from the CHARMM36 (red cross) trajectories in comparison to the corresponding values derived from the experimental R1, R2, and NOE measurements analyzed with global correlation times derived from the bond orientations in the X-ray structure21 (black circles) or from the trajectory-averaged instantaneous effective rotational diffusion rates (blue circles). The assumed experimental uncertainty was 0.5% for R1 and R2 and 0.010 for NOE.

The CHARMM22/CMAP analysis yielded a maximum rmsd value at Gly 14 of only 0.29%, while the maxima for the AMBER ff99SB and ffSB-ILDN analyses were reached for Gly 41 with 0.47% and 0.36%, respectively. This analysis suggested that optimal relaxation analyses may be achieved by setting τN to τM/2 whenever τM is less than 2/ ωN. To test this approach, the autocorrelation calculation and 1 H−15N order parameter analysis was repeated with τM set to 2/ωN. A minimum for the mean and standard deviation of the autocorrelation function fits was observed at this τM value for perturbations on both τH and τN.

rotational diffusion rate for the 1H−15N vectors yielded quite similar Sf2 and SH2 predictions as could be anticipated from their sampling of timeframes significantly faster than global tumbling. On the other hand, the differences between these two approaches for estimating the SN2 values became apparent with a number of the residues. In the X-ray structure, the 1H−15N bond vector for Gly 41 is oriented approximately along the principal axis of the diffusion in the X-ray structure (27.9°) K


Article


experimental 15N relaxation values, the resultant uncertainties estimated for the SN2 values were also 2- to 3-fold larger than for the SH2 values (Figure 10), thus indicating that the Larmor frequency-selective order parameter analysis exhibits a qualitatively similar pattern of responsiveness to both forms of uncertainty in fitting the experimental measurements. As expected from the fact that conventional Lipari−Szabo analysis3 of these data yields an excellent fit to the single parameter S2 model for the majority of backbone amides, most of the backbone exhibits no evidence of internal dynamics in the time frame of either τH or τN. Given the sensitivity of the SN2 values to errors in the effective global correlation time for each vector, as well as the general issue of accuracy in the individual experimental measurements, detailed interpretation of discrepancies for specific bond vectors may at times be problematic. On the other hand, a consistent pattern of differential order parameter predictions provides increased confidence in analyzing these data. As illustrated in Figure 1, the CHARMM36 simulations provide a fairly good prediction of the experimental NMR relaxation values for the modestly flexible loop between the β1 and β2 strands. The similarity within the pattern of order parameters predicted for this segment of Gly 9 through Lys 13 suggests that the experimental dynamics of this loop is somewhat faster than predicted by the CHARMM36 force field (Figure 10). It may be noted that this turn segment is the sole portion of the GB3 backbone for which evidence of conformational dynamics in the approximate μs range has been reported.44 For such conformational transitions which are too fast to yield exchange linebroadening effects but too slow to affect dynamics in the time frame of molecular tumbling, the corresponding metastable states will contribute to the observed relaxation measurements in a simple populationweighted average.

which leads to a comparatively long effective rotational diffusion time. When this highly mobile bond vector is averaged over the CHARMM36 trajectories, the effective internal orientation angle increases to 53.9° which is close to the dipolar magic angle value (54.7°) that is predicted for isotropic sampling. With a rotational anisotropy value of 1.35, this larger effective angle corresponds to a 7.3% shorter diffusion time for this bond vector.32 Since the isotropically averaged rotational diffusion rate is determined by the protein data set as a whole, the order parameter analysis attempts to compensate for this discrepancy by adjusting the SN2 value for the Gly 41 residue. With a more rapidly decaying global autocorrelation function predicted from the CHARMM36derived effective internal orientation angle, the optimized internal dynamics order parameter SN2 is shifted to a considerably higher value (Figure 10). As a result, apparent increased global tumbling dynamics is traded off with apparent decreased internal dynamics, relative to the analysis obtained using the X-ray structure-based bond orientation. The opposite pattern of compensation is illustrated for the dynamics analysis of Asn 37. A quite similar set of observations apply to the corresponding analysis of the experimental GB3 15N relaxation values based on the AMBER f99SB-ILDN force field simulations (Figure S2). The order parameter predictions were initially carried out assuming the isotropically averaged τM value of 3.115 ns which had been obtained by applying the ModelFree algorithm30 to a subset of experimental GB3 data for resonances that were adequately fitted by the single parameter S2 model. Applying this τM value yielded 45 of the SN2 values being set to 1.0, indicative of an artificially shortened global tumbling time. Given the acute sensitivity of the SN2 values to the assumed rotational diffusion rate, the τM values were increased 1.5% and 1.25% (3.162 and 3.154 ns) for the X-ray derived and the CHARMM36-derived bond orientations, respectively. These optimizations were obtained by initially carrying out the order parameter calculations without the constraint SN2 ≤ 1.0. Applying the CHARMM36-derived bond orientations with the optimized average τM value of 3.154 ns yielded apparent SN2 values above 1.0 for the 24 residues. The average margin above 1.0 for these 24 residues was 0.030, which provides an upper bound for the uncertainty in these order parameter estimates. An equal number of residues lay between 0.95 and 1.00, consistent with all of these residues exhibiting no statistically significant internal motion in time frame of τN. A similar analysis was used to optimize the global tumbling time for the analysis of the X-ray derived bond orientations, although a somewhat larger maximum error estimate was obtained (0.033). The apparent uncertainties for the SH2 values of the highly ordered residues could be similarly estimated from the positions with optimized values above 1.0 to yield an uncertainty estimate of 0.010. Excluding seven residues for which either experimental data is lacking or previous MD studies have indicated markedly divergent dynamics predictions, the CHARMM36-derived and experimentally derived values for 48 amides of GB3 agree to an average of 0.016, 0.010, and 0.020 for the fast limit (Sf2) and Larmor frequency-selective (SH2 and SN2) order parameters, respectively (Figure 10). A quite similar level of agreement was observed between the order parameters for this same set of residues predicted from simulation or derived from experimental data for the AMBER ff99SB-ILDN force field analysis (Figure S2). When random variations were applied to the

4. CONCLUSIONS The full exploitation of the synergy offered by joint NMR-MD analysis has arguably been inhibited by the common restriction of simulation studies to the fitting of superimposed trajectory frames to experimentally derived S2 values. In such an approach, not only has the information reflecting the time frame of those dynamics been underutilized, various issues relating to the implicit assumption of separability between global and internal motion have often been neglected. Most notably, for all proteins exhibiting discernible anisotropy in their global tumbling dynamics, the internal bond vector dynamics directly impact both the global and internal autocorrelation functions. We have found that trajectoryaveraged instantaneous rotational diffusion rates markedly improve the accuracy of the resultant global bond vector autocorrelation functions, while rotational velocity rescaling analysis provides a useful upper bond to the errors arising from the separability assumption. For the backbone 1H−15N and the nonmethyl 1H−13C bond vectors of the protein GB3, the wide range of dynamics behavior embodied in the bond vector autocorrelation functions obtained from standard force field simulations are surprisingly well modeled by a simple (Sf2,SH2,SX2) order parameter representation whose parameter values can be directly derived from experimentally accessible R1, R2, and NOE relaxation values. The robustness of this approach to MD-based interpretation of NMR measurements is enhanced by the fact that it can be readily applied in a differential mode which should facilitate the mutual cancellation of residual L


Article


free Approach to the Interpretation of N-15 Nuclear Relaxation of Proteins. J. Am. Chem. Soc. 1990, 112, 4989−4991. (10) Peng, J. W.; Wagner, G. Mapping of Spectral Density-Functions using Heteronuclear NMR Relaxation Measurements. J. Magn. Reson. 1992, 98, 308−332. (11) Ishima, R.; Nagayama, K. Protein Backbone Dynamics Revealed by Quasi Spectral Density-Function Analysis of Amide N-15 Nuclei. Biochemistry 1995, 34, 3162−3171. (12) Farrow, N. A.; Zhang, O.; Szabo, A.; Torchia, D. A.; Kay, L. E. Spectral Density Function Mapping using 15N Relaxation Data Exclusively. J. Biomol. NMR 1995, 6, 153−162. (13) Farrow, N. A.; Zhang, O.; Forman-Kay, J. D.; Kay, L. E. Comparison of the backbone dynamics of a folded and an unfolded SH3 domain existing in equilibrium in aqueous buffer. Biochemistry 1995, 34, 868−878. (14) MacKerell, A. D.; Feig, M.; Brooks, C. L. Improved treatment of the protein backbone in empirical force fields. J. Am. Chem. Soc. 2004, 126, 698−699. (15) MacKerell, A. D.; Feig, M.; Brooks, C. L. Extending the treatment of backbone energetics in protein force fields: Limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics. J. Comput. Chem. 2004, 25, 1400−1415. (16) Huang, J.; MacKerell, A. D. CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data. J. Comput. Chem. 2013, 34, 2135−2145. (17) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. D.; Mittal, J.; Feig, M.; MacKerell, A. D. Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the backbone ϕ,ψ and Side-Chain χ1 and χ2 Dihedral Angles. J. Chem. Theory Comput. 2012, 8, 3257−3273. (18) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins: Struct., Funct., Genet. 2006, 65, 712−725. (19) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the Amber99SB protein force field. Proteins: Struct., Funct., Genet. 2010, 78, 1950−1958. (20) LeMaster, D. M. Larmor Frequency Selective Model Free Analysis of Protein NMR Relaxation. J. Biomol. NMR 1995, 6, 366− 374. (21) Derrick, J. P.; Wigley, D. B. The third IgG-binding domain from streptococcal protein G. An analysis by X-ray crystallography of the structure alone and in a complex with Fab. J. Mol. Biol. 1994, 243, 906−918. (22) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera - A visualization system for exploratory research and analysis. J. Comput. Chem. 2004, 25, 1605−1612. (23) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (24) Anderson, J. S.; LeMaster, D. M. Rotational Velocity Rescaling of Molecular Dynamics Trajectories for Direct Prediction of Protein NMR Relaxation. Biophys. Chem. 2012, 168−169, 28−39. (25) Case, D. A. Calculations of NMR dipolar coupling strengths in model peptides. J. Biomol. NMR 1999, 15, 95−102. (26) Hernández, G.; LeMaster, D. M. Quantifying protein dynamics in the ps-ns time regime by NMR relaxation. J. Biomol. NMR 2016, 66, 163−174. (27) Kearsley, S. K. On the orthogonal transformation used for structural comparisons. Acta Crystallogr., Sect. A: Found. Crystallogr. 1989, 45, 208−210. (28) Rupp, B.; Parkin, S. PDBSUP-A FORTRANP Program That Determines the Rotation Matrix and Translation Vector for Best Fit Superposition of Two PDB Files by Solving the Quaternion Eigenvalue

systematic error. By this means, the dynamical effects resulting from the systematic alteration of force field parameters and their relationship to the corresponding experimental measurements should be more straightforwardly interpretable on a detailed structural basis.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00387. Larmor frequency-selective order parameter fitting to AMBER ff99SB-ILDN H−C autocorrelation functions upon variation of the time parameters (Table S1). Comparison of experimental 15N relaxation values and Larmor frequency-selective order parameters with values derived from AMBER force field simulation (Figures S1 and S2). (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Telephone: (518)474-6396. Fax: (518)473-2895. E-mail: [email protected]. ORCID

David M. LeMaster: 0000-0003-2810-9465 Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS We acknowledge the use of the IBM Intelligent Cluster at Union College. REFERENCES

(1) Levantino, M.; Yorke, B. A.; Monteiro, D. C. F.; Cammarata, M.; Pearson, A. R. Using synchrotrons and XFELs for time-resolved X-ray crystallography and solution scattering experiments on biomolecules. Curr. Opin. Struct. Biol. 2015, 35, 41−48. (2) Mark, P.; Nilsson, L. Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. J. Phys. Chem. A 2001, 105, 9954−9960. (3) Lipari, G.; Szabo, A. Model-Free Approach to the Interpretation of Nuclear Magnetic Resonance Relaxation in Macromolecules. 1. Theory and Range of Validity. J. Am. Chem. Soc. 1982, 104, 4546− 4559. (4) Zeiske, T.; Stafford, K. A.; Friesner, R. A.; Palmer, A. G. Startingstructure dependence of nanosecond timescale intersubstate transitions and reproducibility of NMR-derived order parameters. Proteins: Struct., Funct., Genet. 2013, 81, 499−509. (5) Gu, Y.; Li, D. W.; Bruschweiler, R. NMR order parameter determination from long molecular dynamics trajectories for objective comparison with experiment. J. Chem. Theory Comput. 2014, 10, 2599−2607. (6) Maragakis, P.; Lindorff-Larsen, K.; Eastwood, M. P.; Dror, R. O.; Klepeis, J. L.; Arkin, I. T.; Jensen, M. Ø.; Xu, H.; Trbovic, N.; Friesner, R. A.; Palmer, A. G.; Shaw, D. E. Microsecond molecular dynamics simulation shows effect of slow loop dynamics on backbone amide order parameters of proteins. J. Phys. Chem. B 2008, 112, 6155−6158. (7) Jin, D.; Andrec, M.; Montelione, G. T.; Levy, R. M. Propagation of experimental uncertainties using the Lipari-Szabo model-free analysis of proteins dynamics. J. Biomol. NMR 1998, 12, 471−492. (8) Jin, D.; Figueirido, F.; Montelione, G. T.; Levy, R. M. Impact of the precision in NMR relaxation measurements on the interpretation of protein dynamics. J. Am. Chem. Soc. 1997, 119, 6923−6924. (9) Clore, G. M.; Szabo, A.; Bax, A.; Kay, L. E.; Driscoll, P. C.; Gronenborn, A. M. Deviations from the Simple 2-Parameter ModelM


Article

Journal of Chemical Theory and Computation Problem; Lawrence Livermore National Laboratory, Livermore, CA, 1996. (29) Lakomek, N. A.; Ying, J.; Bax, A. Measurement of 15N relaxation rates in perdeuterated proteins by TROSY-based methods. J. Biomol. NMR 2012, 53, 209−221. (30) Mandel, A. M.; Akke, M.; Palmer, A. G. Dynamics of ribonuclease H: Temperature dependence of motions on multiple time scales. Biochemistry 1996, 35, 16009−16023. (31) Woessner, D. T. Nuclear Spin Relaxation in Ellipsoids Undergoing Rotational Brownian Motion. J. Chem. Phys. 1962, 37, 647−654. (32) Wong, V.; Case, D. A. Evaluating Rotational Diffusion from Protein MD Simulations. J. Phys. Chem. B 2008, 112, 6013−6024. (33) Li, T.; Jing, Q.; Yao, L. Dynamics of the GB3 Loop Regions from MD Simulation: How Much of It Is Real? J. Phys. Chem. B 2011, 115, 3488−3495. (34) Trbovic, N.; Kim, B.; Friesner, R. A.; Palmer, A. G. Structural analysis of protein dynamics by MD simulations and NMR spinrelaxation. Proteins: Struct., Funct., Genet. 2008, 71, 684−694. (35) Lu, C. Y.; vandenBout, D. A. Effect of finite trajectory length on the correlation function analysis of single molecule data. J. Chem. Phys. 2006, 125, 124701. (36) Cavanagh, J.; Fairbrother, W. J.; Palmer, A. G.; Skelton, N. J. Relaxation Mechanisms. In Protein NMR Spectroscopy, Principles & Practice; Academic Press: San Diego, CA, 1996; pp 279, 287. (37) Shaw, D. E.; Maragakis, P.; Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Eastwood, M. P.; Bank, J. A.; Jumper, J. M.; Salmon, J. K.; Shan, Y. B.; Wriggers, W. Atomic-Level Characterization of the Structural Dynamics of Proteins. Science 2010, 330, 341−346. (38) LeMaster, D. M.; Kushlan, D. M. Dynamical Mapping of E. coli Thioredoxin via 13C NMR Relaxation Analysis. J. Am. Chem. Soc. 1996, 118, 9255−9264. (39) LeMaster, D. M. NMR relaxation order parameter analysis of the dynamics of protein sidechains. J. Am. Chem. Soc. 1999, 121, 1726−1742. (40) Khan, S. N.; Charlier, C.; Augustyniak, R.; Salvi, N.; Déjean, V.; Bodenhausen, G.; Lequin, O.; Pelupessy, P.; Ferrage, F. Distribution of Pico- and Nanosecond Motions in Disordered Proteins from Nuclear Spin Relaxation. Biophys. J. 2015, 109, 988−999. (41) Kadeřav́ ek, P.; Zapletal, V.; Rabatinová, A.; Krásný, L.; Sklenár,̌ V.; Ž ídek, L. Spectral density mapping protocols for analysis of molecular motions in disordered proteins. J. Biomol. NMR 2014, 58, 193−207. (42) Hill, R. B.; Bracken, C.; DeGrado, W. F.; Palmer, A. G. Molecular motions and protein folding: Characterization of the backbone dynamics and folding equilibrium of α2D using 13C NMR spin relaxation. J. Am. Chem. Soc. 2000, 122, 11610−11619. (43) Kadeřav́ ek, P.; Zapletal, V.; Fiala, R.; Srb, P.; Padrta, P.; Přecechtĕlová, J. P.; Šoltésová, M.; Kowalewski, J.; Widmalm, G.; Chmelík, J.; Sklenár,̌ V.; Ž ídek, L. Spectral density mapping at multiple magnetic fields suitable for 13C NMR relaxation studies. J. Magn. Reson. 2016, 266, 23−40. (44) Pratihar, S.; Sabo, T. M.; Ban, D.; Fenwick, R. B.; Becker, S.; Salvatella, X.; Griesinger, C.; Lee, D. Kinetics of the Antibody Recognition Site in the Third IgG-Binding Domain of Protein G. Angew. Chem., Int. Ed. 2016, 55, 9567−9570.

N


Prediction of Bond Vector Autocorrelation Functions from Larmor

Recommend Documents