A referencing strategy for the direct comparison of NMR and MD motional parameters in RNA
Nuclear magnetic resonance (NMR) spectroscopy and molecular dynamics (MD) simulations are both techniques that can be used to characterize the structural dynamics of biomolecules and their underlying timescales. Comparison of relaxation parameters obtained through each methodology allows for cross validation of techniques, and for complementarity in the analysis of dynamics. Here we present a combined NMR/MD study of the dynamics of HIV-1 TAR RNA. We compute relaxation constants (R1, R2, and NOE) and model-free parameters (S2 and \xcf\x84) from a 65 ns molecular dynamics (MD) trajectory and compare them with the respective parameters measured in a domain-elongation NMR experiment. Using the elongated domain as the frame of reference for all computed parameters allows for a direct comparison between experiment and simulation. We see good agreement for many parameters and gain further insight into the nature of the local and global dynamics of TAR, which are found to be quite complex, spanning multiple timescales. For the few cases where agreement is poor, comparison of the dynamical parameters provides insight into the limits of each technique. We suggest a frequency-matching procedure that yields an upper bound for the timescale of dynamics to which the NMR relaxation experiment is sensitive.
I. INTRODUCTION<!>A. Relaxation Parameters<!>B. Model-free Parameters<!>C. Experimental Relaxation Parameters<!>D. Model Free Analysis of Experimental Data<!>E. MD simulation<!>F. Calculated Parameters<!>III. RESULTS AND DISCUSSION<!>A. Domain I: The reference domain<!>B. Domain II: Local and collective motions<!>C. Stretched Exponential Behavior<!>D. Discrepancy in the slow motions, and a frequency-matching procedure<!>IV. CONCLUSIONS
<p>In recent years it has become clear that describing not only the structure, but also the dynamics of a biomolecule is essential to its full characterization [1–5]. This is especially important for RNA, for which conformational flexibility is nearly universally observed to be essential for its binding and function [1, 2, 4, 6, 7]. However, obtaining a detailed characterization of the dynamics of RNA can be difficult. NMR spectroscopy is one of the most powerful methods for the detection of internal motions, with atomic resolution and sensitivity anywhere from the picosecond to millisecond timescales for molecules in the solution state [7–9]. However, collecting data on RNA and interpreting dynamical information is not always straightforward, due, in part, to the possible breakdown in assumptions underlying the methods of data interpretation, namely the assumption that the overall and internal motions are decoupled [10, 11]. Molecular dynamics (MD) simulations [12] also provide atomic level detail of solvated molecules and can, in principle, provide a more detailed analysis of structure and dynamics than is possible by NMR. However MD is limited in timescale by the length of the trajectory that can be collected, which can often lead to incomplete sampling of conformational dynamics (i.e., to the so-called broken ergodicity problem [13]). Additionally, there are inaccuracies inherent to the approximations in the empirical force fields employed that can lead to errors in the simulated dynamics.</p><p>The comparison of biomolecular dynamics, as derived experimentally by NMR and theoretically by MD simulations, thus provides an important means for cross-validating and interpreting results from each methodology [58–64, 67–69]. Moreover, these techniques are highly complementary and thus their combined use provides the basis for a synergistic approach to the characterization of dynamics. Comparison of dynamics between MD and NMR is most often made through NMR spin relaxation parameters [14]. This is relatively straightforward when molecular motions are not coupled to overall rotational diffusion. In such a case, the NMR data can readily be interpreted in terms of internal and overall motions [10, 11]. Further, although overall rotational tumbling is usually not sufficiently sampled during the finite trajectory of an MD simulation, it can generally be removed through the overlay of each trajectory snapshot onto a reference structure [15], allowing for the accurate calculation of internal dynamical parameters.</p><p>Comparing results from NMR and MD is significantly more complicated when internal motions are coupled to overall rotational diffusion. Such motional correlations have been detected in some proteins [16–19] and are prevalent in RNA molecules [20–24]. In such cases, parameters describing internal motions are difficult to obtain experimentally because the contributions to relaxation arising from internal motions cannot be easily separated from those due to overall rotational diffusion. Moreover, from a MD simulation standpoint, when these motions widely alter the overall shape of the RNA on the timescale of an MD trajectory, the moment of inertia changes significantly, and the overall structure can no longer be used as a single reference frame for eliminating overall rotational diffusion. Although this problem has been nicely addressed in the iRED analysis presented by Prompers and Bruschweiler, this method does not allow for accurate computation of model free parameters (the most commonly compared NMR parameter) if motions are coupled.</p><p>A recently described domain-elongation NMR experiment addresses the problem of motional coupling specifically for RNA [21]. In this experiment, the coupling of internal and overall motions is largely eliminated through a substantial elongation of one of the helical domains. This experimental strategy slows overall tumbling and causes it to be dominated by the elongated domain, thus decoupling internal motions that were on the timescale of the overall tumbling of the non-elongated RNA, and effectively anchoring the reference frame onto the elongated domain. Here, we mimic this reference frame in the analysis of an MD trajectory by overlaying each trajectory snapshot such that it aligns with the elongated domain, which serves as the reference, allowing us to compute relaxation parameters that can be directly compared to those obtained in a domain-elongation NMR experiment. Using this strategy we examine the dynamics of the wild-type transactivation response (TAR) RNA element of HIV-1 (see Figure 1) for which we have previously shown that the internal and overall motions are inseparable [24]. We find good agreement between many parameters, which reveal a picture of complex dynamics in the RNA spanning multiple timescales.</p><!><p>NMR relaxation parameters report on the overall and internal motions of a molecule. When cross-correlation effects are suppressed, imino 15N relaxation is dominated by dipolar interactions with the directly bonded 1H and by the chemical shift anisotropy (CSA) of the 15N. Assuming an exchange term, Rex, of zero and an axially symmetric CSA tensor with the principal axis co-linear with the NH bond, one can express the longitudinal and transverse relaxation constants (R1 and R2) and the Nuclear Overhauser Enhancement parameter (NOE) as follows [25]: 
(1)R1=d24(3J(ωN)+J(ωH−ωN)+6J(ωH+ωN))+c23J(ωN)
(2)R2=d28(4J(0)+3J(ωN)+J(ωH−ωN)+6J(ωH))+c218(4J(0)+3J(ωN))
(3)NOE=1+d2γH4R1γN(6J(ωH+ωN)−J(ωH−ωN)).</p><p>Here, J(ω) is the spectral density function, ωN and ωH are the Larmor frequencies in rad·s−1, d = (μ0hγHγN/8π2〈r3〉) and c = ΔσωN, with μ0 the permeability of free space, h Planck's constant, γH and γN the gyromagnetic ratios of 1H and 15N respectively, r the NH bond length, and Δσ the CSA of the 15N.</p><p>Relaxation parameters are obtained directly from an NMR experiment and can also be computed from the coordinate trajectory of an MD simulation by calculating spectral density functions from the internal and overall correlation functions. For calculations from a finite molecular trajectory, J(ω) is approximated as</p><p>Here 
Coaxial(t) is the correlation function for overall tumbling when modelled with an axially symmetric diffusion tensor[26], Ci(t) is the correlation function for internal motions, and Ci(∞) is the converged plateau value of Ci(t). The first term in Eq. 4 is computed directly from the MD trajectory with tmax being the maximum time of the calculated correlation functions. The latter term, approximating the unsampled overall tumbling, is calculated from analytical evaluation of the integral.</p><p>When calculating correlation functions from a molecular dynamics simulation of a flexible multi-domain molecule, it is important to define the proper reference frame. Overlaying the atoms of a stable molecular domain across all trajectory snapshots will eliminate all global motions for that domain and will yield correlation functions for bond vectors in that domain which contain information only on the local motions. However, with the same reference frame, bond vectors for another domain in the same molecule will yield correlation functions which contain information on both the local motions as well as motions induced in the bond vector by any inter-domain fluctuations.</p><!><p>The most popular method for interpreting relaxation data in terms of dynamics is the Lipari-Szabo model free approach [27], in which the dynamics of bond vectors are described by amplitudes (order parameter, S2) and corresponding correlation times (τ) as obtained through the parameterization of the correlation function. For simple motions a single exponential fit is used, however for bond vectors undergoing more complex motion, an extended two-exponential form is used [11],</p><p>
(5)Ci(t)=S2+(1−Sf2)e−t/τf+(Sf2−S2)e−t/τs, where 
S2=Sf2Ss2 is the tail value of the time correlation function and the f and s subscripts denote "fast" and "slow" motions respectively. It is assumed that internal and overall motions are not correlated and that the "fast" and "slow" internal motions are also not correlated.</p><p>Model free parameters are obtained from an MD simulation through fitting of the computed internal correlation function, Ci(t). The order parameter, S2, is either obtained as the plateau of the correlation function or is calculated from an equilibrium expression [29] given as</p><p>The corresponding correlation time, τ, is generally obtained from fitting an exponential function to Ci(t) where τ is the time-constant from the fitting. For bond vectors undergoing more complicated motions, a sum of exponentials may be used, providing a series of correlation time constants.</p><!><p>Two samples of elongated wild-type HIV-1 TAR (see Figure 1b) were prepared by in vitro transcription, using a chemically synthesized double stranded DNA template and T7 polymerase. The RNA product was purified on a 15% denaturing polyacrylamide gel and recovered by electroelution and ethanol preciptiation. The RNA was exchanged into pH 6.4 NMR buffer containing 15 mM sodium phosphate, 25 mM sodium chloride and 0.1 mM EDTA. Each TAR sample had 22 Watson-Crick base-pairs added to the terminal end of domain I, for the purposes of obtaining relaxation parameters (Figure 1a). In one sample, domain I was elongated with AU base pairs and transcribed with 13C and 15N labeled guanidine and cytodine. The other sample was elongated with GC base pairs and transcribed with labeled adenine and uradine. This labeling strategy was employed so as to avoid significant peak overlap due to the additional residues in the elongated domain. An overlay of HSQC spectra from the elongated and non-elongated constructs show very few differences, excepting for those residues directly adjacent to the residues added in elongation, indicating that there are very few structural differences between the constructs (see Figure 1d).</p><p>15N relaxation studies were performed on the elongated TAR samples at a field strength of 600 MHz on an Avance Bruker spectrometer equipped with a triple resonance 5 mm cryogenic probe. Longitudinal (R1) and transverse (R2(CPMG)) relaxation rates, and NOEs were measured at 298 K using standard 1D and 2D R1, R2, and NOE experiments on the AU-labeled and GC-labeled samples respectively. In the R1 experiment, the imino 1H spins were decoupled during the relaxation delay using 1250 μs iburp2 pulses centered on the imino signals. A recycling delay of 1.9 s was used, and relaxation delays of 60, 120, 240, 480, 640, 800, and 1200 ms were used, repeating experiments with delays of 120 and 800 ms for the purposes of error analysis. In the R2 experiment a [0013]N phase cycling scheme was employed to supress off-resonance effects[30]. As in the R1 experiments the imino proton spins were decoupled during the relaxation delay using 1250 μs ipurb2 pulses. The inter-pulse delay for the CPMG scheme was 1300 μs. A recycling delay of 1.9 s was used, and relaxation delays of 6.2, 12.4, 24.8, 37.2, 49.6, 62.0, and 74.4 ms were implmeneted, with repeated experiments performed for delays of 12.4 and 62.0 ms for the purposes of error analysis. Rates R1 and R2 were obtained through fitting the relaxation data to exponential functions (with an offset of 0) using the program Origin 7.0 (OriginLab Corp.). The R2 values were corrected for off-resonance effects as described in Yip et. al.[30]. Errors in R1 and R2 were taken to be the fitting errors as deterined using Origin. NOE errors were obtained from NMRDraw and were taken as the square root the sum of the squares of the errors in the two peak heights. Data was obtained for N1H1 or N3H3 bonds in residues G17, G18, G21, and U42 in domain I and residues G28, G36, and U38 in domain II. The peaks for G26 and G43 could not be resolved and thus data for these residues was not obtained. Additionally the lack of a hydrogen bond between A22 and U40 precluded the ability to obtain data for U40.</p><!><p>The relaxation parameters were interpreted according to the model free formalism yielding parameters describing the amplitude (S2) and timescale (τ) of bond vector motion. The program ModelFree provided by the Palmer laboratory (http://www.palmer.hs.columbia.edu/ [31]) was used to best fit the relaxation data to dynamical parameters. The input structure for all calculations was built using ideal A-form models of the two helices and the loop from structure 3 of the NMR solved ensemble of structures of TAR[32], with the average global orientation of the helices obtained from an RDC analysis (data not shown). The program HYDRONMR[33] was utilized to estimate the rotational diffusion tensor parameters based on the hydrodynamic shape of the structure, with input parameters AER = 3.2, NSIG = 6, SIGMIN = 1.4, SIGMAX =2.4, T = 298.15 K, and a solvent viscosity of 0.008906. The resulting diffusion tensor parameters were Dzz = 17.6 MHz, Dyy = 3.53 MHz, and Dxx = 3.47 MHz. The tensor parameters obtained from HYDRONMR were then further refined using ModelFree with the input structure rotated into the determined hydrodynamic principle axis system (PAS).</p><p>For the purposes of the ModelFree calculations the determined errors in R1 and R2 were doubled and the NOE error was set to 8%. A static bond length of 1.01 Å was used [21, 65] and Δσ was set to −130 ppm and −100 ppm for N1 and N3 respectively. All calculations utilized 300 randomly distributed data sets generated by Monte Carlo simulations. To refine the diffusion parameters the rotational correlation time and polar angles θ and φ (defining the orientation of the z-component of the diffusion tensor) were varied from 15–24 ns and from 0–20° and 0–360° respectively, with Dratio set to the HYDRONMR predicted value (Dratio= 5.0). This procedure was performed with data for domain I from the GC-labeled construct utilizing the Powell minimization algorithm, yielding τm = 18.7 ns, θ = 18.1° and φ = 151.2°. These simulations were repeated with 10% variations in Dratio with very similar results. Using a structure in the GC-labeled diffusion tensor frame τm was refined for the AU-labeled domain I data using the simulated annealing protocol with "grid 100 50 0.9", yielding τm= 20.7 ns. Utilizing these parameters the procedure outlined by Mandel et. al.[31] was followed for model selection for individual bond vectors. Model 1 (S2) was chosen for G18 and U42, model 2 (S2 and τ) for G17 and G21, and Model 5 (
Ss2,Sf2, and τs) was chosen for G28, G36, and U38. These parameters were then refined using the Powell minimization and simulated annealing algorithms for the GC-labeled and AU-labeled data respectively.</p><!><p>A 65-ns MD simulation of wild-type HIV-1 TAR (see Figure 1a) was performed using the CHARMM force field with the parameter set 27 [34]. Details of the simulation protocol and the first 20 ns of the resultant trajectory have been reported previously [24, 28, 64]. In what follows, we give a brief description of the set-up for the presently reported trajectory. Starting coordinates were obtained from structure 3 of the family of free TAR NMR structures (PDB 1ANR) [35]. This structure was chosen as it yields the best agreement with previously measured residual dipolar couplings (RDCs) [36]. The RNA was charge-neutralized using sodium counter ions and solvated in a 35 Å sphere of TIP3P water [37], which allowed for > 9 Å distance between the surface of the sphere and all RNA atoms. A stochastic boundary potential on the water molecules was used [38]. This system was minimized and heated to 300 K, while harmonically constraining the heavy atoms of the RNA with a force constant of 62 kcal/mol/Å2 for 100 ps, after which the constraints were removed and the system equilibrated for 1ns. A Nosé-Hoover thermostat [39, 40] was used to maintain a constant temperature of 300 K throughout the simulation, with a 1 fs time step and a bath coupling constant of 50 ps−1.</p><!><p>Relaxation parameters R1, R2, and NOE were calculated from the MD trajectory for N1H1/N3H3 vectors for residues G17, G18, G21, and U42 in domain I and residues G28, G36, and U38 in domain II. In order to be able to rigorously compare the calculated and experimental parameters, it was necessary to match the experimental diffusion tensor frame. Since domain I of TAR was elongated in the experiments, it dominates overall tumbling and is the reference frame for all experimental parameters. Thus, to mimic this reference frame, the MD relaxation parameters were calculated with the domain I heavy atoms as a reference for overlaying trajectory snapshots, excluding A22-U40 as this base-pair is not hydrogen bonded[32]. This referencing strategy anchors the MD reference frame to domain I and should render the computed and experimental parameters comparable. The overall rotational tumbling was modeled by an axially symmetric diffusion tensor, which is supported by the HYDRONMR calculation of the tensor parameters, which found Dxx ~ Dyy. Values for the diffusion tensor were set equal to those derived from the model free analysis of the experimental data, where the principle components of the diffusion tensor, D|| and D⊥, were found to be 19.1 MHz and 3.8 MHz respectively for the Guanine data and 17.3 MHz and 3.5 MHz respectively for the Uracil data, with a magnetic field of 14.1 Tesla, and Δσ of −100 ppm and −130 ppm for N3 and N1 respectively. The domain I axis is, to a good approximation, collinear with the principal axis of the experimental diffusion tensor, thus the average angle of each bond vector with respect to the principal axis of the diffusion tensor was approximated by the average angle between the bond vector and the helix axis of domain I, which was determined using the program CURVES[41]. Values for 〈r3〉 were calculated as an average over the MD trajectory, and were within 0.005 Å of the static bond lengths used in the experimental analysis. All correlation functions for use in calculating the spectral density functions were computed from 1–65 ns of the trajectory with tmax = 6.5 ns (excepting those discussed in section III.E.). For converged correlation functions, Ci(∞) was calculated using the equilibrium expression for S2 (see Eq. 6). However, for those functions which were not converged during the simulation, Ci(∞) was approximated as the tail value of the correlation function, defined as the average value over the last 500 ps in the C(t) dependence.</p><p>Internal correlation functions for domain II bond vectors were calculated using two reference frames. The domain I reference frame yields functions which contain information on both the local motions of the bond vectors within domain II and any motions induced by the fluctuation of domain II with respect to domain I. Additionally, to obtain information on only the local motions within the domain correlation functions were calculated using domain II as the reference frame. Here the terminal residues and those neighboring the bulge and loop were excluded from the overlay, as they are expected to experience significant dynamics. Assuming that the local and domain motions are separable, dividing the correlation functions calculated for bond vectors in domain II when using domain I as the reference (CI(t)) by the correlation functions obtained using domain II as the reference frame (CII(t)) yields correlation functions containing information on only the motions of the bond vector induced by the slow, global motions of domain II with respect to domain I,</p><p>Order parameters describing the local motions of bond vectors were calculated using Eq. 6 (
Seq2) from 1–65 ns of the MD trajectory as well calculated from the tail values 
Stail2 of the correlation functions (see above). To test for the convergence of the correlation functions 
Seq2 values were compared with 
Stail2 values, the functions were considered converged if 
∣Stail2−Seq2∣<0.005. All correlation functions were parameterized by fitting to single, double, and triple exponential functions using non-linear least square fitting of the form</p><p>
(8)Ci(t)=c0+∑j=1ncje−t/τj, where cj, 
j=0,n¯, are constants, n = 1 corresponds to a single exponential, n = 2 to a double exponential, and n = 3 to a triple exponential fit, and where τj is the correlation time constant. The quality of fit was determined by analysis of the χ2/dof and R2 values obtained for each fit[42]. Errors in the correlation functions were approximated as [43]</p><p>
(9)σC(t)=2τint/T(1−C(t)/C(0)) where T is the upper time limit of the trajectory (i.e., T = 64ns) and τint is the integrated correlation time computed as</p><p>Order parameter errors were approximated to be equal to the error in the last point of the correlation function. Errors in the relaxation parameters were estimated from the standard deviation in a set of 1000 computed relaxation parameters in which Gaussian errors were added to each time point of the internal correlation functions with a standard deviation of σC(t). Errors in the correlation time constants were obtained from the fitting errors of correlation functions in which each time point had added to it error with a standard deviation of σC(t). In these fittings the coefficients in the exponential functions were set to those obtained when fitting the functions without any added error.</p><!><p>In an elongated-domain NMR experiment, as described above, one helical domain of an RNA is elongated using a stretch of Watson-Crick base-pairs. This slows down overall rotational diffusion, thus decoupling it from internal motions that occur on timescales approaching the overall tumbling of the non-elongated construct. Additionally, as the elongated domain significantly dominates the rotational diffusion of the RNA it becomes the reference frame for all measured parameters. Thus, dynamical parameters will report only on the local motions of bond vectors in the elongated domain, but for any additional domains will report on both local and inter-domain motions. In this way, elongation can reveal domain motions that are hidden in the non-elongated constructs due to coupling with the overall rotational diffusion.</p><p>A similar frame of reference can easily be obtained in the analysis of an MD trajectory. To accomplish this, the domain which has been elongated in the experiment (though not necessarily in the simulation) is used as the reference frame for overlaying each snapshot from the MD trajectory by minimization of the root mean square atomic displacements of the heavy atoms. As in the experimental approach, this results in any calculated parameters being inherently referenced to this domain, rendering them comparable to the experimental values. We have applied this strategy to HIV-1 TAR RNA (see Figure 1a), for which we obtained relaxation constants R1, R2, and NOE from an elongated-domain experiment and from a 65 ns MD trajectory. Presented below is a discussion of the cross-validation of the relaxation and associated model-free parameters, and an analysis what they reveal about the dynamics of TAR.</p><!><p>Relaxation parameters (R1, R2, and NOE) obtained experimentally and from the 65 ns MD simulation are shown in Figure 2. These parameters are shown as a function of residue and separated by domain. As explained above, dynamical parameters in domain I will report only on local motions, which are expected to be on the picosecond timescale. Since this is significantly shorter than the rate of overall tumbling in the NMR experiment (τm ≈ 20 ns) and than the length of the MD trajectory (T = 65 ns), these motions should be easily resolved by both techniques. Looking at the values for the reference, domain I (residues 17–22 and 40–45, see Figure 2), very good agreement is observed between the calculated (green) and experimental values (gray), with R1 values differing on average by only 3.1%, R2 values by 3.9%, and NOE values by 13.4%. This good agreement between the experimental and computed values in the reference domain shows that the MD simulation is reliably reproducing the internal dynamics in the domain, and also demonstrates the validity of the referencing strategy.</p><p>The local motions in domain I were also characterized through experimental and computed model free parameters. Shown in Figure 3 are the order parameters and correlation times for local motion, again as a function of residue and separated by domain. Excellent agreement between the experimental and computed S2 values is observed for domain I, indicating that the amplitude of bond vector motions gauged by the two techniques is approximately the same, and the average values of 0.88 and 0.85 obtained from NMR and MD respectively signify a level of dynamics expected for a standard A-form helix. [14, 28]</p><p>For the associated timescales the comparison becomes more complicated. Model free analysis of NMR data typically relies on either a single exponential[10] (models 1 and 2 in ModelFree) or, in the extended formalism, a double exponential fit [11] (model 5) of the correlation function, with further parameterization not possible given the limited data available. For all domain I bond vectors either model 1 or 2 was chosen with resulting timescales of <100ps. However, fitting of the MD calculated correlation functions to obtain the correlation times for the motions indicates that more complex motions are present than a single or even double exponential fit can account for. In fact, a triple exponential best fit the correlation functions out of the single, double, and triple exponential fits tested. This was seen previously by Pfeiffer et. al. [44] who also needed to utilize a triple exponential to properly fit their correlation functions calculated from a molecular dynamics simulation trajectory of a protein (and using, importantly, a different force-field for the simulation). The triple exponential fit yields three distinct correlation times. These include an ultra-fast sub-picosecond decay, followed by a fast decay on the ~20–200 ps timescale, followed by an additional slow decay of ~300–7500 ps. The second and third correlation times are shown in Figure 3b for each NH bond. The second correlation time is compared to the experimentally derived constants for the "fast" motions, for which moderate agreement is observed. This fit indicates that there are multi-timescale motions for the NH bond vectors in the simulations, including motional modes that are slower than those observed based on model-free analysis of the experimental relaxation data when assuming a single correlation function. Based on the overall E-TAR correlation time of 20.7 ns, the NMR relaxation data should be sensitive to all three motional modes observed in the MD; however data at a single magnetic field strength may not permit resolution of the multi-timescale nature of the dynamics. Previous studies have shown that low amplitude, long time-scale motions are not reliably detected in model-free analysis [45, 46]. Thus, these timescales are likely simply evading detection in the model free analysis. Overall, the good agreement observed between the calculated and experimental values for these parameters for domain I (see Figure 2) indicates that the timescales sampled in the MD simulation are realistic.</p><!><p>Experimental and computed relaxation parameters (R1, R2 and NOE (experimental only)) for bond vectors in domain II are shown in Figure 2. With domain I as the reference frame, the domain II (Figure 2, residues 26–39) dynamical parameters will report not only on local motions within the domain, but also on any collective motions present with respect to domain I. Computed NOE values for domain II have a high level of error due to the lack of convergence of the correlation functions, and thus are not shown or compared. For both the experimental and computed data sets, the R2 values obtained for the domain II bond vectors are significantly lower than those obtained for domain I, and similarly R1 values in domain II are observed to be significantly higher than domain I. Given the overall tumbling rate of elongated TAR, this pattern indicates that the domain II vectors are experiencing a greater level of dynamics as compared to domain I vectors. However, the difference between the calculated domain I and II values is greater than that observed experimentally, with the domain II R1 values on average 49.1% higher and the R2 values on average 35.2% lower than their experimental counterparts. This indicates that the domain II vectors in the simulation are experiencing a greater level of dynamics than in the experiment, and further that the increased level of domain II bond vector dynamics with respect to domain I bond vectors is even greater in the simulation than is observed experimentally.</p><p>Model free analysis of the experimental data for domain II favored a model with two modes of motion, a "fast" mode and a "slow" mode, for all bond vectors (i.e. model 5). The fast modes have associated parameters similar to the domain I local parameters, excepting G36 which has a lower S2 value. The larger amplitude (S2 = 0.71–0.84) and longer timescales (τ = 1.9–2.1 ns) associated with the slow motional modes, and the similarity across different domain II vectors suggest that perhaps this mode represents motion of the entire domain in a collective manner. However, this cannot be conclusively determined from the NMR data as there is no way to distinguish between local and domain motions in analysis of the experiment.</p><p>In contrast, in the MD analysis, it is possible to separate these motions, as the local and domain dynamical parameters can be computed separately using distinct reference frames, thus allowing for conclusive determination of the nature of the dynamics. Order parameters describing the local motion of bond vectors within domain II were calculated from correlation functions with domain II as the reference frame (CII(t)), i.e., by overlaying each snapshot of the MD trajectory such that the root mean square displacement of only the domain II atoms is minimized relative to a reference snapshot of domain II. These parameters are shown in Figure 3 alongside the experimentally determined parameters for the "fast" motional mode. Importantly, comparing the calculated domain II local order parameters with the calculated domain I order parameters (Figure 3a) shows that, with the exception of G36, very similar amplitudes of motion are observed for both domains. Thus, the increased level of dynamics observed in the domain II bond vectors compared to domain I as detected by the calculated relaxation constants is conclusively not on the local level, but rather collective domain motions must be present to account for the increased level of dynamics. Accordingly, the experimental "fast" and "slow" modes do indeed represent local and domain motions respectively.</p><p>In comparing the experimental and computed S2 values for the local motions, residues G28 and U38 are in very good agreement. However, the computed value for G36 is significantly lower than the experimental counterpart, indicating a larger amplitude of motion. Interestingly, the S2 value computed from only the first 20 ns of the simulation (data not shown) is in much better agreement with the experimental value, indicating that in the latter part of the simulation G36 accesses a more dynamic state. It was recently shown that better agreement was found between experiment and simulation when S2 values were calculated from a set of simulations each started from a distinct structure as accessed through an accelerated molecular dynamics simulation[47], demonstrating that the proper distribution of global conformational states is also necessary for a proper description of local dynamics. Thus, it is quite possible that the discrepancy in the S2 value for G36 is simply due to an improper or incomplete distribution of structural conformations, which is supported by the fact that the correlation function describing the local motions for G36 has not converged over the 65 ns trajectory (see Figure 3c). The conformational state important for the stabilization of G36 is likely related to the cross-loop base-pair between C30 and G34 that has been observed in experiment. When this base-pair is formed, G34 stacks on G36 and stabilizes it to a less dynamic state. This base-pair exists intermittently (~27% of the time) during the first 20 ns of the simulation, but does not form at all in the last 45 ns of the simulation. It is possible that this base-pair is not stabilized sufficiently during the simulation. However, given the length of time during the simulation in which the base-pair closing-opening was sampled, it is also possible that the proper distribution is simply not complete and that the base-pair would reform in a longer simulation.</p><p>As was seen for correlation functions for bond vectors in domain I, a triple exponential was needed to properly fit the correlation functions for local motions in domain II. This yielded correlation times very similar to those observed for domain I, with the exception of G36 for which the slowest decay was almost 3 ns. The second and third correlation times from the fittings are shown in Figure 3b. The second correlation time is compared to the experimental values for the fast mode.</p><p>A second set of correlation functions for domain II was calculated using domain I as the reference frame (CI(t)). These functions contain information on both the local and domain motions. Assuming that these motions are separable, dividing the CI(t) functions by the CII(t) functions yields correlation functions for domain II containing information only on the slow domain motions (Cs(t)) (see Equation 7). These correlation functions, shown in Figure 4c, have clearly not converged to equilibrium during the simulation. However, fitting of these functions can still give some insight into the nature of the slow motions. Order parameters were calculated from the tails of the Cs(t) correlation functions. These values are shown in Figure 4a and are compared to the experimental "slow" S2 values obtained from the model free analysis of the experimental data. The calculated values are 34% lower, on average, than the experimental values, indicating that there is a larger amplitude of domain motion in the simulation than is observed in the experiment. This is consistent with what was observed from the relaxation parameters and demonstrates that the source of the discrepancy between the domain II experimental and computed parameters stems from increased amplitude of dynamics on the global, rather than local, level in the simulation. RDC measurements on elongated TAR also reveal greater amplitude helix motions that observed by spin relaxation measurements likely due to presence of slower motions occurring at nanosecond timescales that can be detected by RDCs but not spin relaxation. Indeed, fitting of the computed functions shows that a triple exponential better describes the function than either a single or double exponential. The second and third correlation times from this fitting are shown in Figure 4b. The third correlation times are in best agreement with the experimental values for the slow mode, however the computed constants indicate slower motions than are observed in the experiment, with the longest computed correlation time 3 times longer than the observed experimental timescale.</p><p>Interestingly, the multi-timescale nature of the global motions observed here through the triple exponential fittings is consistent with what was seen in a previous analysis of the first 20 ns of this simulation [24]. In that analysis, the inter-helical bend and twist were analyzed as a function of time. Both parameters clearly demonstrated periodic fluctuations on both the picosecond and nanosecond timescales. Additionally, it was clear that the nano-second motions were only just beginning to be sampled, which explains the lack of convergence seen in the currently described correlation functions.</p><!><p>In addition to the multi-exponential fitting performed above we also attempted a stretched exponential fitting of the Cs(t) functions. The stretched exponential function is given as</p><p>
(11)C(t)=C(0)exp(−t/τ)β, where 0 < β ≤ 1 and gives a measure of the non-exponentiality. Such a time dependence, also commonly known as the Kohlrausch-Williams-Watts function, is observed frequently in the relaxation of many complex, disordered systems including proteins and other polymers [48]. Here the stretched exponential fittings were statistically comparable to the triple exponential fit and yielded β values ranging from 0.55–0.58 and τ from 10.7–12.7 ns. The stretched exponential fitting indicates that there are dynamics occurring which are best described by an array of correlation time constants. Both an enthalpic and an entropic model have been proposed to explain this disordered behavior. The former is signified by a multi-minimum potential energy surface, whereas the latter arises due to the particular form of the distribution of multiple conformational pathways available around a free energy minimum. A recent molecular dynamics study of NAD(P)H:Flavin oxidoreductase by Karplus and coworkers [49] found that it was a multi-minimum potential energy surface (i.e., an enthalpic model) that produced the stretched exponential behavior for the system studied. Here the dynamics are not sufficiently sampled to concretely draw a conclusion concerning the source of the multi-timescale motions. However, the continuous nature of the dynamics in TAR (i.e., the absence of long-time trapping in energy wells) suggests that it most likely fits an entropic model [50]. Regardless of its origin, the stretched exponential behavior along with the lack of convergence of the Cs(t) correlation functions lead us to conclude that there are even longer timescales and more complex domain motions in TAR not yet sampled in the current simulation.</p><!><p>For the wild-type TAR system studied here, the discrepancy between the experimentally detected slow motions and those observed in the MD simulation could be due to several factors. One possibility is that the simulation is simply modeling too much dynamics. This could be due to either the nucleic acid force field or the solvent force field. The reliability of the nucleic acid force field is supported by the accurate representation of the internal motions in domains I and II, though this does not conclusively prove its overall accuracy. Concerning the solvent force field, it is a well established fact that none of the current models accurately reproduce all of the properties of bulk water. For the TIP3P model, which was used in this study, discrepancies are observed in the kinetic properties, with a viscosity that is ~ 30% lower than that observed experimentally and a diffusion constant ~ 240% higher [51]. This has been shown to affect the overall rotational diffusion of molecules and it has been suggested that it could also affect the dynamics of solvent exposed residues [52], and thus, it may have an effect on inter-domain motions. The overall rotational diffusion of molecules can be adjusted by a simple constant, however the effects of the viscosity on intramolecular motions is likely not so straightforward as it would be determined not only by interaction with the solvent but also by the biomolecular dihedral angular transitions and other intrachain interactions ("internal friction"). However, the equilibrium properties of the TIP3P water model accurately reproduce experimental properties, thus it is unlikely that this would affect the equilibrium amplitude of inter-domain motions, but rather only the correlation time constant, thus this cannot completely explain the observed discrepancy.</p><p>Additionally, it has been shown that the collective helix motions are strongly dependent on ionic strength [53], decreasing as the ionic strength increases. The level of sodium in the simulation approximates what was in the buffer used in the NMR relaxation experiments, however it is possible that incomplete modeling of the ion/RNA interaction in the simulation could also lead to some discrepancies in the inter-domain motions.</p><p>A more likely possibility is that there are larger amplitude motions with timescales approaching the overall rate of molecular tumbling of the elongated construct that are not fully captured by the NMR relaxation data, but that are being accessed in the MD simulation. If there were motions going undetected it would be similar to the spectral density function containing only a portion of the correlation function for the bond vector in question. To simulate the affects of such a scenario, we re-calculating the relaxation constants from the MD simulation using correlation functions with increasingly truncated values of tmax. The best agreement between the domain II experimental and computed R1 and R2 values was obtained for tmax= 1.8 ns (~ 0.1 × τm). These recomputed relaxation parameters are shown in Figure 5 with the experimental parameters. Compared to the values calculated for tmax = 6.5 ns, the recomputed domain II R1 and R2 values are in much better agreement with the experimental values, differing on average by only 15%. This reduced timescale sensitivity may arise because the fast collective internal motions may re-set the timescale sensitivity limit for other motional modes, i.e. motions that are slower than a dominant fast motional mode but slower than overall tumbling may not significantly affect measured relaxation data, with the added complication that the modes may not be fully decoupled from one another [24]. Importantly, only negligible changes are observed in the domain I values using tmax= 1.8 ns vs. 6.5 ns. Additionally, the S2 values calculated from the Cs(t) correlation functions with tmax= 1.8 ns are in very good agreement with the experimentally observed "slow" S2 values (see Figure 5d), and it is also noteworthy that this timescale (i.e., 1.8 ns) is in very good agreement with the experimentally observed timescale of the slow motions, which is observed to be 1.9–2.1 ns.</p><p>We also computed S2 values for C-H bonds throughout TAR and compared them with values previously reported by Hansen et. al. [66] for UUCG-TAR, see Figure 6. The values shown were computed from CI(t) correlation functions with tmax=6.5ns (Figure 6a) and tmax=1.8ns (Figure 6b). Very good correlation is seen between the experimental and computed values for the lower part of domain I. However, for the more flexible residues (i.e. the upper part of domain I, the bulge and domain II) poor correlation is seen with tmax=6.5ns. However, decreasing tmax to 1.8ns leads to a much better correlation between the experimental and computed values.</p><p>Altogether the MD data strongly suggests that there are larger amplitude, longer timescale motions in TAR. Previous RDC analysis (sensitive to motions up to the millisecond timescale) of UUCG-TAR indeed gives evidence for significantly larger amplitude motions [36]. In that analysis, fitting the RDC data to a cone motional model yielded an amplitude of motion (cone angle) of ~45°, which would correspond to a "slow" S2 value of 0.35. This is a much larger range of motion than what was found through the relaxation analysis, which provides an average slow S2 of 0.76, corresponding to a cone angle of ~24°. Indeed this is an even larger range of motion than in is being accessed in the MD analysis presented here, which yields a "slow" S2 of 0.5, suggesting further that there are even longer timescales and more complex domain motions in TAR not yet sampled in the current simulation. Longer trajectories may allow for a quantitative description of these motions, though likely alternative approaches based on enhanced dynamical sampling methods [54–57] will be necessary to access relaxation kinetics on these long time-scale motions.</p><!><p>The problem of finding a well defined single reference frame has led to great difficulty in calculating dynamical parameters that can be compared to experimental counterparts for flexible molecules such as RNA. However, an elongated-domain NMR experiment on RNA leads to a readily defined reference frame for dynamical parameters, and by mimicking this reference frame in the analysis of an MD trajectory we were able to quantitatively compare MD and NMR relaxation data for wild-type TAR RNA. In this study we obtained relaxation and model-free parameters from an elongated domain experiment and from a 65 ns simulation. Our results reveal that TAR undergoes complex dynamics on the local and global levels, spanning multiple timescales. Further, the MD data indicate the presence of slower, larger amplitude motions than are observed in the NMR relaxation studies, which is in agreement with preliminary RDC analysis that strongly suggests there are indeed larger scale collective dynamics that are evading detection by the relaxation techniques. This self-consistent referencing strategy provides a more rigorous comparison of dynamical parameters than is afforded by traditional techniques when dealing with highly flexible systems such as TAR, and should facilitate the synergistic development of MD and NMR for the characterization of RNA functional dynamics.</p>