it is now widely recognized 
 that many rna molecules are predisposed 
 to forming complexes with proteins by fluctuating spontaneously through 
 an ensemble of structural states .
protein 
 recognition is referred to as  conformational capture  . 
 a description of the physical principles involved in forming rna  protein 
 complexes via conformational capture requires complete description 
 of the dynamics of these structurally labile rna molecules , including 
 a characterization of long- and short - lived conformational states 
 sampled by the rna . however , experimental characterization of transitory 
 states is complicated by the fact that the rate of transition may 
 be too fast to allow for a comprehensive catalogue of all states . thus , although progress has been made recently 
 in experimentally isolating the partially folded states of proteins , the complete elucidation of protein or nucleic acid dynamics requires 
 analytical and computational modeling to complement experimental observations . 

a common approach toward this goal is to fit a limited set of model 
 parameters to experimental data that are sensitive to dynamics , such 
 as nmr relaxation rates , line shapes , or residual dipolar couplings .
this procedure involves guessing - and - checking , using physical constraints 
 on possible motions of the labeled residue(s ) to guide the model - building 
 process .
however , this semianalytic approach becomes complicated when 
 further degrees of freedom and new free parameters are required by 
 the model to fit the data adequately .
a purely computational 
 approach to the description of molecular 
 states requires an accurate potential energy function ( pef ) , followed 
 by either molecular dynamics ( md ) or energy - minimization calculations . 

molecular dynamics simulations are able to generate dynamic trajectories 
 of the molecule and , in principle , explore molecular parameter space 
 if sufficient numbers of trajectories are available ; examples are 
 found using the amber and charmm packages and others , such as lindorff - larsen et al .
however , the extrapolation of dynamics in nucleic 
 acids to time scales of the order of microseconds or longer , where 
 many conformational changes are expected to take place , has only recently 
 begun to be explored .
energy - minimization techniques 
 also rely on a well - validated energy 
 function but involve the subsequent alteration of the relative conformations 
 of parts of the molecule in an iterative manner to find the global 
 energy minima .
it is then possible to generate multiple sample structures more 
 easily than for a complete md calculation . in the current manuscript
, 
 we utilize structures generated by energy - minimization techniques 
 as a complement to md - based analyses of dynamics .
use of energy - minimized 
 structures is facilitated by the availability of structures on the 
 rosetta server .
other servers such as 
 the mc - sym / mc - fold pipeline are also 
 available that allow the user to obtain all - atom rna models .
we use the hiv-1 tar ( trans - activation response ) rna molecule as 
 a model for a dynamics simulation based on a set of 500 low energy 
 models generated using the program farfar for the 29-nucleotide apical section of the rna ( figure 1 ) .
the tar rna binds the 
 viral regulatory protein tat , a critical 
 transcription elongation factor essential for viral replication .
the 
 tat binding site surrounds the single - stranded , trinucleotide ( ucu ) 
 bulge and is contained within the 29-nucleotide 
 construct .
the bulge region interlinking the two helical stems is 
 the primary binding site for the tat protein and exhibits considerable flexibility allowing for the two helical regions to adopt a wide range of relative 
 orientations .
, where the free tar 
 rna exchanges between multiple conformers , one or more of which are 
 amenable to tat binding . tar rna also 
 provides a common rna structural motif , where two helices are connected 
 by a single - stranded bulge on one end and a backbone  hinge  
 at the other . other rna s exhibit similar structures , such 
 as k - turns , or the hiv-1 rre ( rev - response 
 element ) .
it is therefore worthwhile to characterize 
 the dynamics of such fundamental motifs , especially to characterize 
  large - scale  motions ( as opposed to more localized 
 motions ) of one helical domain relative to another .
29-nucleotide apical 
 stem - loop of hiv-1 tar rna : ( a ) secondary 
 structure ; ( b ) sample tertiary structure .
the different submotifs 
 and the component residues used in the analyses are color - coded : upper 
 helix ( turquoise ) , lower helix ( red ) , and the single - stranded bulge 
 ( green ) .
earlier work in our group used 
 solid - state nmr to identify intermediate 
 rate motions for the residue u38 in the upper helix .
we more recently characterized the motions of u40 and u42 
 in the lower helix ( wei huang , unpublished data ) .
these data indicate 
 that the two helices move relative to each other at slow rates ( 10 to 10 s ) relative to the 
 rotational diffusion rate of the entire molecule ( 10 s ) .
these results motivated us to look at the 
 distributions of the orientations of the upper helix relative to the 
 lower helix among the set of lowest energy structures , as characterized 
 by the euler angles of an upper helix - attached reference frame relative 
 to a different frame attached to the lower helix . reduction to a set 
 of three angles characterizes many features of the dynamical trajectory 
 of this particular structural motif , while simplifying a very large - dimensional 
 problem to a set of three essential coordinates .
local base librations , 
 i.e. , rotations of the bases around a base normal ( representing vibrations 
 of the base around the equilibrium base - paired orientation ) or rotations 
 of the base around the glycosidic bond ( for the single stranded bases ) , 
 are included in the simulations as well and provide atomistic detail 
 regarding the motions of individual residues . here
we extend 
 our prior studies of nonrigid rotation and helix 
 reorientation in tar rna to a full description of the concerted dynamics 
 of the upper and lower helices and the bulged loop .
our approach describes 
 a dynamic trajectory based on an ensemble of energy minimized structures 
 with the rosetta program farfar . because farfar does not provide a 
 boltzmann - weighted distribution of states ( and
thus does not provide 
 the entropy and free energy ) our calculations rely on fitting phenomenological 
 parameters to the data . to reduce the complexity of the problem , a 
 predominant set of conformers
we used experimental residual dipolar couplings 
 ( rdcs ) obtained from partial alignment of the tar rna by one alignment 
 medium to achieve this selection , and other groups have recently used 
 similar filtering techniques .
in addition to filtering out a set of long - lived 
 conformers along the stochastic trajectory , we acquire information 
 regarding their relative populations .
we complement the rdc filter 
 with principal component analyses ( pcas ) based on a judiciously chosen 
 set of backbone torsion angles , as done to establish conformational 
 clustering and free energy landscapes of rnas and proteins , and a pca can serve as a useful aid in setting up dynamics calculations .
finally , 
 we use the outcome of that analysis to calculate solution 
 spin  lattice relaxation times t1 and the rotating - frame 
 spin  lattice relaxation time t1 for c spins in nucleotides located in the upper and lower helices 
 and in the bulge .
the protocol would allow extensions to other observables 
 as well , but the current work focuses on these to provide a clear 
 proof of method .
we direct the interested reader to more extensive 
 work on other nmr relaxation parameters and their simulations .
the first approach 
 ( slow exchange or  se  model ) assumes an effectively 
 infinite time scale of exchange ( i.e. , the conformational exchange 
 is much slower than any other relevant motional time scale ) . the second 
 approach ( general rate or 
gr  model ) allows for an 
 arbitrary time scale of exchange ( including one that overlaps with 
 the rotational diffusion time scale ) . by combining complementary experimental 
 and analytical techniques into a single framework
, we have been able 
 to construct a viable dynamic trajectory for the tar rna .
protocol 
 used in this manuscript for the simulation of domain - dynamics 
 in a molecule .
we present here the solid - state nmr ( ssnmr)-derived models 
 that served as motivations for the solution relaxation simulations . 

ssnmr studies of the dynamics of the uridine bases u38 , u23 , and u25 
 in tar were carried out using samples selectively deuterated at the 
 5- and 6-carbon base sites .
u38 was chosen to represent 
 the dynamics of the upper helical stem , whereas u23 and u25 are of 
 interest on account of their positions in the single - stranded bulge . 

more recently , 5,6-h labels were introduced at the u40 
 and u42 positions in tar , in correspondence with the lower helix ( u42 ) 
 and an unstable base pair that closes the bulge region ( u40 ) ( wei 
 huang , unpublished data ) .
the data included line shapes as well as 
 t1z and t1q relaxation times collected on samples 
 hydrated in all cases to 16 water molecules per nucleotide to reproduce 
 conditions where motions were shown to be solution - like .
each of the upper helical and single - stranded sites investigated 
 had characteristically distinct spectral features . within the lower 
 helix ,
u40 and u42 showed results similar to each other , but distinct 
 from the upper helix or single - stranded sites .
motional models generated 
 to fit the data included a slower base motion consisting of jump between 
 two equally populated sites , superposed on a faster motion occurring 
 around the normal to the base - plane for the helical residues and around 
 the glycosidic bond for the bulge residues .
the analyses of the spectra 
 recorded for the upper and lower helical sites were done independently 
 resulting in two different sets of parameters .
the u38 base 
 was modeled as undergoing a two - site 4 jump process around 
 the base - normal at a rate of 2  10 s in addition to a two - equivalent - site conformational exchange process , 
 where the upper helix underwent a 9 bend and a 15 twist 
 at a rate of 10 s relative 
 to a crystal - fixed frame .
a similar model 
 could fit the data for the lower helix bases u40 and u42 as well , 
 resulting in a 018 bending motion and a 1825 
 twisting motion of the lower helix , but at slower rates on the order 
 of 10 s. by analyzing relaxation data , 
 we observed small amplitude ( 6 to 9 ) local 
 motions of the base at a rate of 10 s for both u40 and u42 .
the u23 data were fit by significantly different 
 models , involving a local two - site jump about the glycosidic bond 
 of 11 at a rate of 10 s , and a 24 hop of the base at a rate of 10 s , whereas u25 was modeled as experiencing a 
 30 jump at 6  10 s in 
 addition to a much slower twisting of the base ( 6  10 s ) with a large amplitude of 40. we utilized the same two - site jump models for the local base motions 
 of both the helical and bulge residues , under the rationale that solid - state 
 sample conditions should be able to replicate solution conditions 
 at least in the small - amplitude local motions at the hydration conditions 
 of our studies .
we also used the time scale of the local base jumps 
 as a starting point for simulations of the solution conditions .
finally , 
 we used the observation that conformational exchange motions in solution 
 conditions occur on a time scale similar to that in solid - state samples . 

formalism for part of the relaxation simulations , which assumes an 
 exchange process much slower than all other motional rates .
we did , 
 however , also consider the more general case of an arbitrary rate 
 of conformational exchange .
the torsion angles that were altered in the generation of the structures 
 used here were those of the residues in the bulge and the bulge - adjacent 
 base pairs of the two helical domains .
the 500 lowest energy structures 
 represented a distribution in energy of about 13 rosetta units , where 
 1 rosetta unit is approximately 1 kbt ( rhiju das , private communication ) . to reduce the 
 dimensionality of the problem , we assumed that the simple motif of 
 a helix  bulge  helix can be represented by a reduced 
 set of relative helical orientation parameters .
a set of three euler 
 angles transforming between the two helices is taken to be sufficient 
 to discriminate between helical conformations ( figure 3 ) if the helices behave as rigid units
. relation 
 between the euler angles determined from the atomic coordinates 
 and the associated domain motions . we have used experimental residual dipolar couplings ( rdcs ) 
 for 
 several bonds along the molecule . because rdcs are long - time weighted 
 averages over all conformational states of the molecule , they potentially 
 provide means of extracting both the best - fit conformers and their 
 relative populations .
the structures that were selected by this rdc - filtering 
 process , along with the best - fit populations , were subsequently used 
 to simulate the solution relaxation times .
we incorporated the 
 local base motions obtained by fitting the 
 ssnmr data into the rdc - selected structures and calculated the solution 
 relaxation times using previously published methods .
when a small amount of variation for the local amplitudes and rates 
 of motion , and for the rates of conformational exchange between the 
 selected structures , was allowed , the solution relaxation times for 
 the 17 residues in the molecule could be calculated ( 7 in the lower 
 helix , 3 in the bulge , and 7 in the upper helix ) with this method .
the relative 
 orientations of the upper and lower helices were quantified 
 by considering the orientation of the normal vector of the u38 base 
 relative to a frame defined by a z - axis aligned with 
 the lower helical axis .
the method of evaluating this relative orientation 
 and the subsequent binning of structures within this scheme is as 
 follows:(a)define 
 the upper and lower helical 
 axes for all structures using the program 3dna .
the upper helical axis is taken to be the average local 
 helical axis of the a27-u38::g28-c37 and g28-c37::c29-g36 dinucleotide 
 steps , where the base pairs flanking the bulge were excluded due to 
 possible distortions from ideal a - form helical structure .
the lower 
 helical axis is calculated similarly as the average over the c19-g43::a20-u42 
 and a20-u42::g21-c41 dinucleotide steps.(b)define the lower helix coordinate 
 frame ( lhf ) by choosing the lower helical axis ( calculated above ) 
 as the z - axis and the perpendicular from the z - axis to the g43 c8 atom as the y - axis . 

this choice of the y - axis was made as the 500 structures 
 did not differ in the orientations of the first few base pairs ( including 
 c19-g43 ) , so the y - axis would be the same across 
 all structures .
however , the specific choice of the g43 c8 atom for 
 this purpose was arbitrary.(c)calculate the  angles for each 
 structure , defined as the angle between the projection of the upper 
 helical axis onto the xy - plane of the lhf and the x - axis of the lhf.(d)calculate the  angles , defined 
 as the angle between the upper helical axis and lower helical axis.(e)the  angle is then 
 defined 
 by the orientation of the normal of u38 base ( the vector perpendicular 
 to the c4c2 and c6c2 bonds ) about the upper helical 
 axis . extracting this information from the structures
requires first 
 removing the  and  dependence by rotating the original 
 u38 base - normal vector vu38norm about the fixed lhf axes as follows : vu38norm = rylhf()rzlhf()vu38norm .
the resultant vectors are distributed 
 around the lhf z - axis as a function of their  
 angles .
the euler angles described above are related to the domain 
 motions as shown in figure 1.(f)bin the 500 structures as a function 
 of the euler angle set {  ,  , }. the bins were 
 chosen in 10 increments for  ( 0   
  360 ) and  ( 0 
  
 180 ) angles , resulting in 36 and 18 bins , respectively . instead 
 of binning the  angle in terms of degree increments , we fixed 
 the number of bins to 5 because of an observed correlation between 
 the  and  angles , which results in a shift in the values 
 of  for every  bin , as well as the restriction of the 
  values to only a portion of the full phase space .
thus , trying 
 to bin all possible  angles would have unnecessarily increased 
 the computational time .
define 
 the upper and lower helical 
 axes for all structures using the program 3dna .
the upper helical axis is taken to be the average local 
 helical axis of the a27-u38::g28-c37 and g28-c37::c29-g36 dinucleotide 
 steps , where the base pairs flanking the bulge were excluded due to 
 possible distortions from ideal a - form helical structure .
the lower 
 helical axis is calculated similarly as the average over the c19-g43::a20-u42 
 and a20-u42::g21-c41 dinucleotide steps .
define the lower helix coordinate 
 frame ( lhf ) by choosing the lower helical axis ( calculated above ) 
 as the z - axis and the perpendicular from the z - axis to the g43 c8 atom as the y - axis . 

this choice of the y - axis was made as the 500 structures 
 did not differ in the orientations of the first few base pairs ( including 
 c19-g43 ) , so the y - axis would be the same across 
 all structures .
however , the specific choice of the g43 c8 atom for 
 this purpose was arbitrary .
calculate the  angles for each 
 structure , defined as the angle between the projection of the upper 
 helical axis onto the xy - plane of the lhf and the x - axis of the lhf .
calculate the  angles , defined 
 as the angle between the upper helical axis and lower helical axis .
the  angle is then 
 defined 
 by the orientation of the normal of u38 base ( the vector perpendicular 
 to the c4c2 and c6c2 bonds ) about the upper helical 
 axis . extracting this information from the structures
requires first 
 removing the  and  dependence by rotating the original 
 u38 base - normal vector vu38norm about the fixed lhf axes as follows : vu38norm = rylhf()rzlhf()vu38norm .
the resultant vectors are distributed 
 around the lhf z - axis as a function of their  
 angles .
the euler angles described above are related to the domain 
 motions as shown in figure 1 .
bin the 500 structures as a function 
 of the euler angle set {  ,  , }. the bins were 
 chosen in 10 increments for  ( 0   
  360 ) and  ( 0 
  
 180 ) angles , resulting in 36 and 18 bins , respectively . instead 
 of binning the  angle in terms of degree increments , we fixed 
 the number of bins to 5 because of an observed correlation between 
 the  and  angles , which results in a shift in the values 
 of  for every  bin , as well as the restriction of the 
  values to only a portion of the full phase space .
thus , trying 
 to bin all possible  angles would have unnecessarily increased 
 the computational time . to select a subset 
 of structures for multisite jump simulations from a starting set of 
 500 energy - minimized structures generated by farfar
we only consider the data 
 from a peg / hexanol mixture because this 
 uncharged medium aligns the charged rna only through steric hindrance 
 of the overall rotation .
this situation does not require a detailed 
 characterization of the rna charge density , as would be required for 
 the simulation of rdcs measured in charged alignment media .
extension 
 to charged media data may be considered in the future . the purely 
 steric version of pales was used to calculate the eigenvalues and 
 eigenvectors of the saupe alignment tensor of the molecule by sampling 
 multiple allowed orientations in the presence of a flat ( infinitely 
 large ) obstruction .
the alignment tensor eigenvalues are then used 
 in conjunction with directional information for a particular bond 
 relative to the alignment tensor principal axis frame ( pas ) . for the 
 purposes of simulating the experimental rdcs ,
the only motions 
 to be considered are then the exchanges between distinct conformers . 
 to include dynamic sampling of multiple conformations , we use the 
 weighted average of the rdcs of each conformer rdcistotal = i=1nconformerspirdcisi for the rdc of a bond between spins i and s , with 
 the population of each conformer i being represented by pi .
the residues considered for 
 this study were c19 , a20 , g21 , a22 , 
 u40 , c41 , u42 , and g43 from the lower helical stem , u23 , c24 , and 
 u25 from the bulge , and g26 , a27 , g28 , c29 , g36 , c37 , u38 , and c39 
 from the upper helical stem .
the bond types included the c6h6 
 bond ( pyrimidines ) , the c8h8 bond ( purines ) , and the c5c6 
 bond ( pyrimidines ) from the bases , the c1n1 ( pyrimidines ) 
 and the c1n9 ( purines ) glycosidic bonds , the c1h1 
 and c4h4 bonds from the furanose rings , and 
 the c5h5 and c5h5 
 bonds from the backbone .
furthermore , the rdcs for the bulge residues 
 were considered separately , and not in the same simulations as the 
 helical residues .
this is because the bulge is likely to be significantly 
 more mobile than the helices and will sample a significantly larger 
 number of configurations , requiring a different set of simulation 
 conditions .
effective bond lengths of the species in question 
 are required 
 in the final pales calculation .
the values we chose for the bond lengths 
 are the aliphatic c  h bond length , rch , aliph = 1.1  , the aromatic c  h bond 
 length , rch , arom = 1.09  ( average 
 of ying et al . and allen et al . ) , the aromatic c  c bond length , rcc , arom = 1.4  , the glycosidic c  n bond length for the cytosines , rcn , glyccyt = 1.47  , and the glycosidic c  n bond lengths for the remaining 
 base types , rcn , glycrest = 1.48  . to gain a geometric perspective on the 
 best - fit structures from 
 the rdc comparison , we binned the simulated couplings for each structure 
 according to the {  ,  ,  } angles determined for 
 the 500 structures .
the rdcs within each bin were then uniformly averaged 
 to produce a single bin - rdc for each bond type and residue , rather 
 than retaining the values for individual structures .
this was done 
 because , although the structures within each bin have similar helical 
 conformations , they may differ in the orientations of bonds of certain 
 residues relative to the large - scale conformations . by averaging over 
 these variations in the bonds
, we effectively included in the simulations 
 the small amplitude thermal fluctuations of the atomic bond orientations . to select the best - fit set of structures , we started with an arbitrary 
 initial set of n structures and allowed the choice 
 of the ( n + 1)th angular bin to float while attempting 
 to optimize the total  as well as the pearson s 
 correlation coefficient between the simulated and experimental rdcs . 

in addition , we varied the relative weights of the ( n + 1 ) bins . once the best - fit parameters were obtained for this first 
 iteration ,
the ( n + 1)th bin so obtained replaced 
 the lowest probability structure from the remaining n , and the search was repeated for a second iteration and beyond . 

to make the final choice of n , we started with n = 2 and after minimizing the  for a 
 choice of two structures , we added a third and repeated the process . 

proceeding thus , we found that n = 5 structures ( i.e. , 
 5 bins which also happened to have only one structure in each of them ) 
 were sufficient to produce the best - fit to the rdc data , with no improvement 
 after the fifth iteration .
in addition , as a separate , independent 
 check of the results of the previous procedure and to allow for the 
 variation of the number of bins more easily , we generated a markov - chain 
 monte carlo ( mcmc ) simulation to search through the bins and populations 
 for a best - fit to the rdcs , with the potential for varying the number 
 of bins .
the number of parameters for an n - conformer 
 search was 2n  1 ( n bin 
 choices + n  1 probabilities to be floated ) . 

the markov - chain monte carlo method did not yield better results than 
 the iterative technique . to determine the number of exchanging conformers 
 required to describe 
 the dynamics in tar rna , and to corroborate the results of our rdc 
 filter
, we performed a principal component analysis ( pca ) , following 
 procedures applied to molecular ensembles of proteins , and rnas . the covariance matrix ij = ( qi  qi)(qj  qj ) 
 was calculated as described by mu et al . ,
i.e. , by proposing the following variable set { q2j}:1where j is the jth torsion angle of interest , 
 and ntorsion is the number of torsion 
 angles used in the analysis
. the use of the cosine and sine functions 
 removes complications associated with the periodicity of the torsion 
 angles by helping to uniquely identify particular values of the angles .
the covariance matrix of these variables is calculated with averages 
 over the full ensemble of structures and is subsequently diagonalized . 

the eigenvalues ( and associated eigenvectors ) are arranged in descending 
 order , with the highest values representing modes with the largest 
 contributions to the structural scatter . in our results , we have found 
 that 2 or 3 modes contain a majority ( 70% ) of the total variance 
 in the data .
in addition , histograms of the projections of the ensembles 
 for each eigenvalue are examined for gaussianity . as discussed in 
 the above references , a mode whose histogram consists 
 of a single gaussian - like peak only represents continuous fluctuations 
 about a central structure , whereas multimodal distributions describe 
 discrete conformations separated by free - energy barriers .
thus , from 
 the perspective of assessing the conformational transitions of the 
 molecule , only modes with multimodal distributions are considered 
 relevant .
our results for pcas with different choices of torsion angle 
 sets ( as described in the following ) indicate that the first 2 or 
 3 modes were non - unimodal in distribution and therefore of primary 
 significance in describing conformational exchanges .
the first set included all torsional  suites  along the bulge and hinge region ( pca method 
 1 ) .
the sugar - to - sugar suite as defined in the above reference consists 
 of the set of 6 backbone torsion angles {  ,  ,  , 
  ,  ,  } as well as the glycosidic torsion angle 
 . in our study , we focused attention only on the 6 backbone 
 torsion angles .
the residues included were a22 , u23 , c24 , u25 , g26 , 
 c39 , and u40 .
these were meant to encompass the conformationally relevant 
 part of the molecule under the assumption of relatively rigid helices . 

the number of torsion angles considered was therefore 42 ( 7 bases 
 6 torsion angles ) resulting in an 84 dimensional pca ( because 
 the cosine and sine functions of the angles are treated as separate 
 variables ) .
however , the pca of this set resulted in several modes 
 which contributed substantially , with no clear clusters in any of 
 the largest modes .
the inclusion of the single - stranded bulge 
 region may have caused 
 the lack of structure in the pca projections calculated in this first 
 method , because single - stranded regions are significantly more flexible 
 and may add a level of disorder to the torsional distribution . to 
 isolate the interhelical motions , we evaluated the pca of only those 
 torsional  suites  that extend from u38 to c39 and from 
 c39 to u40 ( pca method 2 ) .
these suites include only the hinge region 
 of the molecule , and resulted in a 24-dimensional pca .
to check 
 the robustness of the clustering results obtained from 
 method 2 , we then carried out a series of pcas by successively including 
 an additional sugar - to - sugar backbone suite : u38 through c41 ( pca 
 method 3 ) , u38 through u42 ( pca method 4 ) , c37 through u42 ( pca method 
 5 ) , and c37 through u40 ( pca method 6 ) .
the results of these analyses 
 and the relation to the results from the rdc fit will be discussed 
 in the results below .
the evaluation 
 of the solution relaxation times is based on two previously explored 
 techniques .
both methods 
 involve the calculation of the two - time correlation function for the 
 orientation of an atomic bond located within a nonrigid brownian rotator 
 relative to a fixed laboratory frame .
the evaluation proceeds by introducing 
 an intervening reference frame associated with the principal axis 
 system ( pas ) for the rotational diffusion tensor , which is time dependent 
 as a result of exchange between different structural conformers .
the 
 correlation function then becomes dependent on wigner rotation matrices 
 whose arguments are the euler angles that orient the rotational diffusion 
 tensor of the molecule relative to the laboratory frame . in turn , 
 the fokker 
planck equation in three - dimensions allows the evaluation 
 of the transition probabilities from one set of euler angle orientations 
 to another .
planck equation , which accounts 
 for coupling between rotational diffusion and conformational changes , 
 is2where p(,,t|0, ) is the probability 
 that the molecule will transition from a diffusion tensor - to - laboratory 
 frame orientation of 0 at time 0 and in 
 a conformational state  to an orientation of  
 and a conformational state  at time t , given 
 diffusion tensor elements of dij in the 
 conformational state  and an exchange rate of r between conformational 
 states  and . this equation can also formally cover 
 the case of continuous transitions to new conformational states by 
 allowing the number of conformers to be infinite .
however , in all 
 cases considered in the original references , and here as well , only 
 discrete jumps will be considered .
if the exchange rate is much 
 slower than the rate of diffusion and all other time scales of exchange 
 processes , yet faster than the relaxation rates themselves , the rotational 
 diffusion problem is effectively decoupled from conformational exchange . 
 under these conditions ,
the se model applies , and it is possible to 
 calculate the relaxation rates for the molecule as the weighted - average 
 over the relaxation rates for the individual conformational states 
 of the molecule:3 expressions for the 
 noes can also be derived by appropriate weighting 
 of the noes of individual conformational states .
the procedure to 
 calculate the relaxation times for each conformer are provided in 
 the original reference .
solution of eq 2 for the case of a general 
 exchange rate ( i.e. , the gr model ) between conformers involves considerably 
 more analytical and computational processing .
we 
 have solved this problem for the case in which the eigenvectors of 
 the diffusion tensors of the exchanging conformers are coincident 
 at the moment of exchange ( see ryabov et al . for general case ) and correlation functions have been published .
this general rate analysis also indicates that 
 the slow exchange regime occurs at time scales longer than about 1 
 s . for the carbons considered in this analysis ,
the rotating 
 frame z - relaxation time t1 was 
 measured instead of t2 . under 
 the application of a weak spin - lock field and the assumption of lorentzian 
 spectral densities
therefore , for the purposes of this 
 paper we operate under this assumption and simulate the t2 relaxation times . the structure set and the corresponding 
 populations preselected 
 by the rdc filter are the basis for simulations of the relaxation 
 times using eq 2 .
eigenvalues and eigenvectors 
 of the rotational diffusion tensor are first calculated using the 
 public - domain program hydronmr
. then 
 orientations of the atomic bonds of the residues of interest are calculated 
 with respect to this axis system , i.e. , the principal axis system 
 of the rotational diffusion tensor .
the orientational parameters , 
 together with the diffusion tensor eigenvalues are input into the 
 two algorithms we have derived for simulations of the relaxation times 
 of nonrigidly rotating macromolecules : ( a ) the  slow - exchange  
 formalism , describing the case where the 
 conformational jumps occur at a rate much slower than the rate of 
 overall rotational diffusion of the molecule , and ( b ) the general 
 rate formalism , where arbitrary rates of exchange 
 are allowed .
the  slow - exchange  formalism , though merely 
 a limiting case of the general rate theory , has the advantage of being 
 significantly faster and easier to implement and so is considered 
 here .
the residues considered were a20 , g21 , a22 , u40 , c41 , 
 u42 , and 
 g43 from the lower helical stem , u23 , c24 , and u25 from the bulge , 
 and g26 , a27 , g28 , c29 , g36 , u38 , and c39 from the upper helical stem . 
 in the current work
, we have only simulated the motions of the bases 
 of these residues : the c6h6 bonds for the pyrimidines and 
 the c8h8 bonds for the purines .
the parameters used 
 in the simulation include the atomic element 
 radius ( aer ) , the radius of the beads used 
 in the hydronmr calculation of the diffusion tensor , and the bond 
 lengths .
the aer was chosen to be 2.3  and the bond lengths 
 for the carbon 
hydrogen bonds for the aromatic bases were chosen 
 to be 1.1  , both choices having been justified in previously 
 published work . the viscosity was chosen 
 to be 1.096 cp to correspond to the conditions 
 of the solution experiments ( 99.9% d2o at 25 c ) .
we have also incorporated the two - site base motions inspired by 
 simulations of the solid - state nmr ( ssnmr ) data : the so - called  base 
 libration  occurs around a vector normal to the plane of the 
 base in the case of helical residues , whereas the two - site motion 
 is modeled to be around the glycosidic bond for the bulge residues . 

we floated the values of the rates and amplitudes of these two - site 
 jumps relative to the ssnmr models , which were found to be on a time 
 scale much shorter than that of the conformational exchange .
for the 
 slow exchange simulations , these internal , local motional rates and 
 amplitudes were the only free parameters , whereas in the general rate 
 simulations we floated the conformational exchange rates between the 
 states .
the fitting procedures were carried out in a combination of 
 grid - searches and mcmc techniques .
five structures 
 provided the best  values to fit 
 the rdc data , as obtained by iteratively searching through the bins 
 and updating the choices of bins and relative populations .
they are 
 shown in figure 4 , and key characteristics 
 are summarized in table 1 . for simplicity , 
 the structures will be referred to by their respective bend angles . 

thus , the highest population structure will be called the  45 
 structure  , the second highest populate structure as the 
five structures obtained by the rdc - filtering 
 procedure to represent 
 the ensemble of tar conformation that describes the experimental data , 
 together with their interhelical euler angles and population percentages .
the  for the best - fit 
 set of structures was 
 11  460 for a set of 48 rdcs , with 9 degrees of freedom ( 5 bin 
 choices and 4 probabilities ) , for a reduced r2 = 302 ( = 11460/(48 
  9  1 ) ) , whereas the pearson s correlation 
 coefficient was 0.72 .
it is to be noted that the  represent unweighted values ; i.e. , the discrepancies between the 
 experimental and simulation rdc values are not inversely weighted 
 with the error bars of the rdc values ( which have not been calculated ) . 

this assumption is equivalent to assuming that the error on all measurements 
 is 1 hz , which is likely to be a considerable underestimation of the 
 true error bars .
the comparison between experiment and simulation 
 is shown in figure 5 . in the figure ,
we further attempted an mcmc fit procedure with the possibility of 
 6 , 7 , or 8 structures but were unable to improve upon the fit .
it 
 is possible that this fit may be improved by the continuation of the 
 mcmc procedure to a greater numbers of iterations , but the pca analysis 
 reported in the following section provides further corroboration that 
 the model found by rdc fitting represents the conformational landscape 
 of the molecule well .
comparison of the experimental rdcs ( red triangles ) with 
 the rdcs 
 generated by the best - fit simulation parameters ( blue circles ) for 
 the helical residues in hiv-1 tar rna .
the values shown include those 
 for backbone ( c5c5 , c5h5 ) , 
 furanose ( c1h1 , c4h4 ) , 
 glycosidic ( c1n1 for pyrimidines , c1n9 
 for purines ) , and base ( c5c6 and c6h6 for pyrimidines , 
 c8h8 for purines ) bonds .
plotting the calculated rdcs against the experimental values 
 ( figure 6 ) , we found that the trend , on average , 
 is toward 
 an underestimation of the rdcs by the simulations .
the dashed blue 
 line in figure 6 is the best - fit line to the 
 data and has a slope of 0.4 and a y - intercept of 
 3.9 .
an underestimation of the rdcs may arise from using a 
 smaller degree of alignment in the simulations than in the actual 
 experimental situation .
one possible source of this discrepancy may 
 be the current assumption of a simple steric model for the alignment 
 of the molecule by peg / hexanol .
recent work has proposed that there are subtleties in the alignment process , 
 including the possible contributions of complex alignment medium topology 
 and electrostatic alignment , that are not incorporated in simulations 
 using only the basic steric version of the pales algorithm
. plot of calculated 
 rdc values vs experimental rdc values for the 
 best - fit set of structures and populations .
the dashed blue line is 
 a linear fit to the data , and the solid red line is the ideal case 
 where the calculated and experimental values match perfectly .
this result is a cautionary statement 
 in the application of simple 
 steric models in the simulation of potentially complex alignment media . 

we attempted to fit rdc data collected in glucopone / hexanol mixtures 
 and in pf1 filamentous bacteriophage media and found that the models 
 selected were different ( two or three of the structures chosen were 
 the same compared to the peg / hexanol model ) .
one obvious reason for 
 this was the availability of rdc data for different bonds in the different 
 media .
however , there is potentially a fundamental difference in the 
 alignment properties of the media as well .
for example , the pf1 phage 
 medium is negatively charged and , thus , 
 as may be the case with peg / hexanol as well , the alignment has an 
 electrostatic component .
this must be taken into account more carefully 
 in future analyses . to further discern the causes of the discrepancy 
 in the fit
, we 
 looked at all helical rdc values that contributed to a deviation of 
 magnitude 10 or greater and found a set of 18 rdcs : 8 correspond to 
 c1h1 furanose ring bonds , 2 to c4h4 
 furanose ring bonds , and the remaining 6 to c6h6/c8h8 base bonds , 
 occurring in 12 helical residues in both the upper and lower stems 
 as well as among all four nucleotide types .
the large deviations observed 
 for furanose ring bonds indicate the existence of additional motions 
 localized to the furanose rings , as has been reported for dna , that have not been accounted for using the current 
 set of 500 structures .
these motions involve an exchange between the 
 c2-endo and c3-endo conformations , and on the time 
 scale of the rdcs these motions may be averaged out to produce an 
 intermediate conformation .
it is likely that the furanose ring samples 
 these conformations even for the residues stacked in a helical configuration . 

a similar argument holds for those base bonds that show a large discrepancy 
 in rdc values : there may be additional vibrations in the base orientations 
 that are not adequately sampled by the 500 structures .
we also 
 compared the five structures selected above to the set 
 of 20 lowest energy tar rna structures recently generated on the basis 
 of noe data that have not been constrained 
 by rdc data . upon calculating helical axes and orientations in the 
 same manner as for the 500 farfar structures in this manuscript ,
we 
 find that 11 out of 20 of the structures , including 3 out of the 5 
 structures that best fit the noe data , occur within the same 10 
 bend angle bin as the highest population structure from the rdc fit 
 described above .
thus , we believe that our approach identifies a predominant 
 conformation set . to test the robustness of the search algorithm , 
 we performed two 
 simulations of fitting a reduced data set upon the random removal 
 of ( a ) 10 rdcs and ( b ) 15 rdcs ( different rdcs were deleted in each 
 of these two cases ) . removing
the first set of 10 resulted in the 
 selection of the same 5 conformers as from the full set , along with 
 a sixth new structure with a population of 4% .
the populations of 
 the 5 full - set best - fit structures were slightly different ( maximum 
 change of 12% ) .
removal of 15 rdcs reproduced 4 out of the 5 full - set 
 best - fit structures , and two new conformers with populations of 7% 
 and 5% .
the maximum change in population among the 4 best - fit structures 
 was 7% in this case .
these results indicate that the choice of structures 
 is robust to a reduction in the size of the experimental data set .
jumps between the five conformers shown in figure 4 require a combination of bending and twisting about either 
 the lower or the upper helices . among the entire set of energy - minimized 
 structures
, there was an observed correlation between the  
 and  angles , as has been reported previously .
these correlations may be reflected 
 even in the jumps among this set of five structures , representing 
 a free energy landscape where the exchanges between minima involve 
 coupled shifts in euler angle values .
we carried out a series 
 of pcas to ( a ) identify the choice of torsion angles that best captures 
 interhelical motions , ( b ) corroborate the rdc - filtered set by overlaying 
 the five chosen structures on the clusters obtained from the pca of 
 choice , and ( c ) identify jump matrix elements for the exchanges between 
 the five rdc - filtered structures in the dynamics calculation .
given the lack of clear clustering results from pca method 1 , which 
 incorporates all backbone torsion angles in the bulge and hinge , we 
 examined whether torsion angles on one side of the helical joint would 
 suffice to describe the interhelical reorientation ( pca method 2 , 
 incorporating the backbone torsion angle suites between u38 and u40 ) . 

when this was done , only the first two principal components ( pcs ) 
 contributed significantly , accounting for 75% of the fluctuations 
 in the molecule ( figure 7 ) .
furthermore , these 
 two pcs were the only ones with a multimodal probability distribution 
 across the 500 structures .
this non - single - gaussian distribution signals 
 the presence of conformational clusters in energy minima separated 
 by significant free energy barriers .
the map of 
 these two pcs is shown in figure 7a and shows 
 the presence of three to four major conformational clusters .
the large 
 red dots superimposed on the 2d plot correspond to the five structures 
 selected from the rdc filtering procedure .
parts b and c of figure 7 show the population distribution histograms for 
 pc 1 and pc 2 , respectively .
again , the positions of the five structures 
 from the rdc filter are indicated by red arrows .
principal component analysis 
 of the torsion angle  suites  
 along the hinge ( residues u38 through u40 ) .
( a ) principal components 
 1 and 2 for the 500 structures , represented by blue dots , with the 
 positions of the five best - fit structures marked explicitly by red 
 dots , and further encircled for clarity .
( b ) histogram of the first 
 principal component ( corresponding to the largest eigenvalue of the 
 covariance matrix ) , with the positions of the five best - fit structures 
 marked explicitly by red arrows .
( c ) histogram of the second principal 
 component ( corresponding to the next - to - largest eigenvalue of the 
 covariance matrix ) , with the positions of the five best - fit structures 
 marked explicitly by red arrows .
first , the five structures chosen 
 as the best - fit set match up well with the main conformational clusters 
 obtained from this pca , suggesting that our rdc - filtered set has captured 
 the relevant information about the major conformational clusters of 
 the molecule .
second , structures with similar interhelical bend angles 
 have similar values of each of the principal components .
the principal 
 component values of the two structures with bend angles of 115 
 and 132 occur in close proximity to each other . the same is 
 true of the pair of structures with bend angles of 61 and 76 
 ( it is true , however , that the 45 structure does not differ 
 significantly from the 61 structure ; the pca suggests that there 
 is a free energy barrier that separates even these two neighboring 
 structures ) .
this is important because other pca methods ( described 
 below ) separate structures with similar interhelical orientations , 
 possibly due to the presence of additional degrees of freedom that 
 do not contribute significantly to interhelical reorientation .
finally , 
 the sums of the probabilities of the best - fit structures within each 
 cluster are similar to each other : the 45 bend structure has 
 a population of 29% , the 61 and 76 structures have a joint 
 population of 40% , and the 115 and 132 structures have 
 a joint population of 31% . because the histogram heights do not correlate 
 well with the nearly uniform probability distribution , we fit the 
 jump rates numerically , as described below .
subsequent attempts 
 at testing the robustness of our cluster analysis 
 yielded further interesting results .
pca method 3 ( u38 through c41 ) 
 showed a marked difference in clustering of the structures .
three 
 pcs contributed significantly , yielding about 70% of the total fluctuations , 
 with all three now being multimodal .
however , the five structures 
 do not all seem to fall within major clusters .
moreover , the 61 
 and 76 structures no longer fall within the same cluster , nor 
 do the structures with bend angles of 115 and 132. this 
 occurs because of the intervention of the torsion angles associated 
 with u40 .
when the helical parameters of the 500 structures were searched , 
 about 37% of the structures u40 did not form a canonical watson 
crick 
 pairing with a22 , indicating considerable conformational variability 
 of this residue , at least within the physical picture generated by 
 the energy - minimization ensemble . to test this hypothesis , we extended 
 the pca up to u42 ( pca method 4 ) on the lower stem and up to c37 on 
 the upper stem ( pca method 5 ) .
both these methods gave very similar 
 clustering to method 3 , with three significant pcs .
we interpret this 
 observation as confirmation that , after the inclusion of u40 , the 
 remaining residues behave fairly rigidly and do not change the results 
 of the pca . as a final confirmation of this conclusion
, we carried 
 out a pca including only torsion angles from c37 to u40 ( pca method 
 6 ) .
the clustering results for this pca proved to be the same as for 
 pca method 2 . given the change in clustering associated with 
 the inclusion of 
 the u40 degrees of freedom
, we utilized the results of pca method 
 2 to set up the rate matrix for the relaxation time simulations . in 
 general , for n conformers , the number of combinations of pairwise 
 rate constants that need to be fit is c2 .
given the clustering suggested by 
 the pca , we reduce the fit problem from 10 ( = c2 for the five rdc filter structures ) to five parameters 
 in the following manner . for pairs of structures that occur within 
 the same cluster in the two pc distributions , we allow for only one 
 distinct exchange rate between both members of the pair in that cluster 
 and any structure in another cluster .
thus , the rates used 
 in the fitting process were for the following exchanges : these exchange processes are 
 shown graphically in figure 8 .
the assumption 
 in this parameter reduction is 
 that all the structures within a cluster are separated from other 
 clusters by similar free energy barriers .
exchange between clusters 
 inferred from the pca , together with 
 the five rates used in the jump matrix for the relaxation time simulations .
the rdc - filtered conformer set 
 of five structures , together with 
 the relative probabilities , were used to calculate the t1 , t2 and noe values .
we used 
 two different approaches to calculate the c6h6 ( pyrimidine ) 
 and c8h8 ( purine ) relaxation times : ( a ) the slow exchange 
 method , where the assumption is that the exchanges occur at an infinitely 
 slow rate ( compared to the rotational correlation time ) , and ( b ) the 
 general rate method , where we fit the relaxation times by allowing 
 the exchange rates to vary arbitrarily .
the general rate 
 analysis has shown that the slow exchange 
 regime in tar rna effectively occurs for time scales longer than about 
 1 s .
thus , the results of this subsection assume conformational 
 exchanges occur on a scale longer than 1 s . as a starting point 
 for the fitting process , we use rates and jump amplitudes for the 
 base librations close to those obtained to fit the solid - state nmr 
 data of the uridines , and 
 changed both the rates and amplitudes in small increments to improve 
 on the . using the fit to the t1 values 
 as a benchmark , we found that it suffices to fit only two base - libration 
 rates , one for the upper helix residues and one for the lower helix 
 residues .
this simulation model was inspired by the results of the 
 two solid - state nmr analyses .
the exceptions , however , are the parameters 
 for u40 and u42 in the lower helix .
we obtain a significantly better 
 fit by using the rates and jump amplitudes for the upper helix for 
 these two residues , indicating that these residues are more similar 
 in local base motion to the upper helix .
 windows  
 in the  plot for the local motion parameters , with 
 the upper and lower helical parameters behaving independently .
the 
 rates in these  windows  vary between 10 s and 10 s , whereas the amplitudes are less than 20. we did attempt 
 simulations of additional models of base - libration such as a treatment 
 of the rates of purines and pyrimidines independently , and the assumption 
 of a constant rate across the entire molecule .
there is always the 
 possibility that more complex models of base libration rates may fit 
 the data better ; for example , we may treat the libration rate of each 
 individual residue as an independent parameter , or choose the purines 
 and pyrimidines in the lower helix as independent from the purines 
 and pyrimidines , respectively .
however , this would increase the number 
 of free parameters in the problem , and we chose the above model as 
 a balance between an arbitrary increase in free parameters and an 
 attempt at a physically realistic representation .
the following 
 representative values of the local motion parameters 
 simulate the relaxation times well:(a)upper helix , u40 and u42 base - libration 
 rate = 4.6  10 s(b)upper helix , u40 and u42 base - libration 
 jump amplitude = 13.7(c)lower helix ( without u40 and u42 ) 
 base - libration rate = 6.6  10 s(d)lower helix ( without 
 u40 and u42 ) 
 base - libration jump amplitude = 9.8 upper helix , u40 and u42 base - libration 
 rate = 4.6  10 s upper helix , u40 and u42 base - libration 
 jump amplitude = 13.7 lower helix ( without u40 and u42 ) 
 base - libration rate = 6.6  10 s lower helix ( without 
 u40 and u42 ) 
 base - libration jump amplitude = 9.8 the match between the experimental and simulation relaxation t1 and t2 values 
 is shown in figure 9 . for quantitative comparison , 
 we calculated the root - mean - square deviation ( rmsd ) across the 14 
 helical t1 values to obtain an rmsd of 
 5.3 ms .
relaxation time simulations for the c6h6 ( pyrimidine ) and 
 c8h8 ( purine ) bonds using the slow exchange method , and comparisons 
 of residuals ( discrepancies relative to experimental values ) to statistical 
 error bars : ( a ) t1 simulations ( blue circles ) compared 
 to the experimental t1 values ( red triangles ) ; ( b ) t2 simulations ( blue circles ) compared to the experimental t2 values ( red triangles ) ; ( c ) difference between simulation 
 t1 and experimental t1 values , together with 
 the statistical error bars on experimental data ( red dashed lines ) 
 at 3.2 ms ; and ( c ) difference between simulation t2 and experimental t2 values , together with the statistical 
 error bars on experimental data ( red dashed lines ) at 0.5 ms .
the error bars shown in figure 9 describe 
 the statistical error in the measurements .
however , the potential 
 systematic error , not quantified by the authors , is larger . yet , even at this level of uncertainty there 
 is a consistent difference between the best - fit local jump amplitudes 
 for the upper and lower helices , with the lower helix amplitudes being 
 larger on average ( this is not reflected in the parameter set shown 
 above but is true for the best - fit windows in general ) .
thus , 
 the solution relaxation times can be fit to an rmsd close 
 to the statistical error in the experiments by assuming a slow exchange 
 rate ( i.e. , slower than 10 s ) between 
 the five structural configurations of tar rna shown in figure 4 with populations determined by rdc filtering and 
 further corroborated by principal component analysis .
the model assumed 
 base - libration rates that range between 10 and 10 s , and libration - rate - dependent jump 
 amplitudes less than about 20. the solid - state nmr parameters 
 for the local base - libration parameters fall within the windows described 
 above , consistent with the notion that solid - state experiments are 
 able to capture the solution - state local motions accurately .
although the slow 
 exchange method provides an approximate scale for exchange processes , 
 the general rate method could further resolve the values of the exchange 
 rates between conformers .
for ease of comparison , we used the local 
 base motion rates and jump amplitudes from the slow - exchange fit , 
 whose results were shown in figure 9 .
although 
 several combinations of 5 pca - inspired rates were found that could 
 be fit to the 14 helical t1 and 14 helicalt2 values with comparable rmsd s , in 
 all cases the inter- and intracluster exchange rates were on the order 
 of 1010 s , which 
 confirms the validity of the slow exchange approximation .
for example , 
 the set of inter- and intracluster exchange rates in table 2 yield the base relaxation times shown in figure 10 with an rmsd for t1 values of 5.2 ms and an rmsd for t2 values 
 of 1.4 ms .
interestingly , there is no clear distinction between the 
 intercluster and intracluster rates , as may be expected from a significant 
 difference in the free energy barriers .
relaxation time simulations for the c6h6 
 ( pyrimidine ) and 
 c8h8 ( purine ) bonds using the general rate method with inter- 
 and intracluster exchange rates from table 2 , and using the same local base motion parameters as for figure 9 : ( a ) t1 simulations 
 ( blue circles ) compared to the experimental t1 values ( red triangles ) ; ( b ) t2 simulations ( blue circles ) compared to the experimental t2 values ( red triangles ) ; ( c ) difference between 
 simulation t1 and experimental t1 values , together with the statistical error 
 bars on experimental data ( red dashed lines ) at 3.2 ms ; ( c ) 
 difference between simulated and experimental t2 values , together with the statistical error bars on experimental 
 data ( red dashed lines ) at 0.5 ms . to summarize ,
if we fix the local motion parameters to the 
 best - fit 
 set obtained using the slow - exchange method but float the conformational 
 exchange rate parameters , we obtain a slight improvement in the quality 
 of the fit .
best fits to experimental relaxations using general rate 
 theory are achieved with inter- and intracluster exchange rates on 
 the order of 1010 s , thus justifying the slow exchange approximation of the prior section . 

we did not float both local motion parameters and conformational exchange 
 parameters simultaneously , but it is reasonable to assume that the 
 results will be similar , especially when we constrain the scale of 
 the local base motion parameters to those observed under solid - state 
 conditions .
in addition to helical 
 residues , we also simulated the relaxation times for the bulge residues 
 using the slow exchange and general rate algorithms . in table 3 , we present relaxation time simulations assuming 
 exchange between the conformations shown in figure 4 using the slow exchange and general rate algorithms .
the 
 base rotation rate is selected as 5  10 s and the amplitude as 15 ( to yield good matches to the t1 values ) .
for the general rate simulations , 
 inter- and intracluster exchange rates shown in table 2 were assumed .
we obtain simulated t1 values that 
 are within 2% of experimental data , but simulated and experimental t2 values for u23 deviate by about 16% ; for c24 , 
 the relative deviation is even greater .
this is likely due to the 
 fact that the t2 relaxation time has a 
 spectral density component ( the j(0 ) term ) that makes 
 this observable sensitive to motions much slower than the larmor frequencies 
 of carbon-13 and protons ( which are on the order of nanoseconds ) . 

the fact that we have been unable to capture the t2 values may indicate that there are additional slower 
 motions of these relatively underconstrained residues that are missing 
 from the conformer set we have selected .
although we have been 
 able to successfully match most of the solution 
 relaxation times to almost within the statistical error bars using 
 the rdc - filtered conformer set , we must address the basic question , 
 are the relaxation times sensitive to motions occurring at rates on 
 the order of microseconds or slower ? though the t1 time has no spectral density dependence slower than 
 a time scale of a nanosecond , the t2 are 
 determined by slower motions as well and their expressions contain 
 a dependence on the j(0 ) spectral density . more importantly , 
 the slow exchange and general rate methods depend on the fact that 
 the time scales of rotational diffusion of many molecules ( including 
 the tar rna considered here ) overlap with the time scales to which 
 both t1 and t2 are sensitive .
thus , two different conformers , with slightly different 
 diffusion tensors will have different characteristic relaxation times 
 when calculated separately .
even the slow exchange averaging process 
 will result in a unique linear combination of relaxation rates that 
 becomes discernible when enough data points are compared .
the general 
 rate theory loses some sensitivity to dynamics for rates much slower 
 than a microsecond , but our least - squares simulations for relaxation 
 times have shown that there is still discernible information to be 
 gained at these time scales .
in this manuscript , we introduce a methodology based on energy - minimized 
 structures that ties together structural and dynamic data , as well 
 as solid - state and solution nmr , to build a dynamic trajectory for 
 the hiv-1 tar rna to atomic - level detail .
the protocol uses ( a ) solid - state 
 nmr data to acquire information about local motions of the bases , 
 ( b ) solution nmr rdc data to identify conformational states and their 
 relative populations , ( c ) pca analyses to identify degrees of freedom 
 relevant to the overall reconfiguration of the molecule , as well as 
 to corroborate the clustering of structures and choose parameters 
 for the dynamics analysis , and ( d ) solution nmr relaxation time simulation 
 techniques previously developed to simulate experimental data and 
 fit the jump rates between the molecular conformers .
the experimental 
 data utilized covers a wide spectrum of motional time scales , from 
 the picosecond scale for solution relaxation times and micro- to nanosecond 
 times derived from solid - state nmr line shapes to the submillisecond 
 time scale investigated with rdcs .
the method reproduces relaxation 
 data at multiple helical residues within the molecule using only five 
 structures out of a set of just 500 possible conformers .
the 
 advantage of this approach lies in the coverage of multiple 
 time scales , including long time scales that are difficult to sample 
 with md , and the ease with which energy - minimized conformers may be 
 obtained for small - to - midsized molecules using public - access software 
 like the rosetta suite of programs .
the 
 use of these structures inverts the problem of dynamics relative to 
 md methods : instead of starting with an initial structure and running 
 the clock forward from t = 0 , we filter out long - lived 
 structures using experimental data and interconnect them in a stochastic 
 trajectory .
thus , the results from such a technique would prove valuable 
 to corroborate md - based simulations , and the protocol could provide 
 a less computationally intensive alternative to extended molecular 
 dynamic simulations .
we have been able to simulate relaxation 
 times of most of the helical 
 residues in the molecule with limited conformational sampling from 
 an already small set of structures .
the success in matching the experimental 
 data indicates that , to gain an understanding of the gross motional 
 properties of rna , it is sufficient to sample the limited phase space 
 of the particular structural motifs that constitute the molecule . 

this protocol was designed to derive a geometric picture of interhelical 
 reorientation based on the limited conformational space available 
 to the torsion angles in and around the bulge and hinge regions .
the 
 assumption that the helices are rigid and the requirement to close 
 the loop formed by the single stranded bulge , adjacent helical base - pairs 
 and the hinge backbone should restrict conformational possibilities 
 for the entire molecule .
steric hindrances and limitations on the 
  stretchability  of the single - stranded region would 
 further impose constraints on molecular reorientation . in practice
, we found that parameter 
 space for the reorientation of one helix relative to the other was 
 expanded by the possibility of one of the bulge adjacent base pairs 
 opening up . among the energy - minimized structures , a significant number 
 had a missing a22-u40 base pair
the u40 base often forms a stabilizing 
 interaction with u25 instead and , among the full set of 500 structures , 
 sometimes with u23 and c24 as well .
we allowed these possibilities 
 to occur in our sample set to reflect fluctuations in the residue 
 orientations , as well as the impacts of these fluctuations on overall 
 molecular conformation .
thus , we believe that models generated using 
 well - vetted potential energy functions can identify sites where new 
 intramolecular bonding and conformational variability might occur .
furthermore , we made a conscious choice to characterize the structural 
 bins by their euler angles .
this choice of parameters has been made 
 previously to enhance reproducibility and 
 comparability to other analyses of the molecule .
such a parametrization 
 represents the core ingredient of most molecular analyses : reduction 
 of the dimensions of the problem to render it tractable .
molecular 
 studies often aim to distill out a few degrees of freedom that are 
 implicated in determining either the structure or dynamics of the 
 system , and many different techniques ( ramachandran plots , phenomenological models , pca 
 analyses ) are directed toward identifying a minimal 
 set of relevant coordinates .
the methodology 
 relies on the assumption that energy - minimized structures sufficiently 
 populate the available conformational space , i.e. , on the assumption 
 of ergodicity .
if ergodic behavior holds , then a sufficiently representative 
 characterization of the energy landscape of the phase space will allow 
 a calculation of the requisite time averages of observables . in the 
 case of the tar rna , the structural motif ( helix  bulge  helix ) 
 is fairly simple , and it is possible to cover a large region of the 
 interhelical orientation space with a relatively small number of structures . 
 for more complex structural motifs , it would be necessary to generate 
 a sample set that covers both the space of molecular reorientations 
 and the range of conformations of local residues relative to a fixed 
 large - scale molecular orientation .
even in the current work , it is 
 possible that we have under - sampled the full range of conformations 
 available to the bulged loop .
a richer sample of both interhelical 
 orientations and orientations of bulge residues relative to particular 
 interhelical orientations may improve the rdc fit .
a second 
 assumption is that the slow exchange ( se ) and general rate ( gr ) methods 
 assume coincidence of the diffusion tensors at the moment of exchange . 

this is not a significant problem for conformers that are not significantly 
 different , but it could pose problems for conformers that are widely 
 separated in conformational space .
we have not attempted to quantify 
 the actual deviation in relaxation times on account of this assumption . 

finally , it bears mentioning that there is a 2-fold 
 degeneracy 
 in the choice of the unit eigenvectors of the rotational diffusion 
 tensors , with the negative of a given choice of unit vector being 
 acceptable as the eigenvector as well .
the choice of eigenvectors 
 does not change the results for the slow exchange formalism ( the expressions 
 are invariant to such changes ) but does impact the general rate theory 
 expressions .
for example , keeping the z - axis the 
 same , the two choices of a right - handed coordinate system vary by 
 180 and would artificially introduce such an extra jump into 
 the calculations .
the means of consistently dealing with such a jump 
 is to track the diffusion tensors as a function of changes in the 
 shape of the molecule , either visually or geometrically and ensure 
 that there is no additional change in the diffusion tensor orientations 
 due to axis inversions . in our particular case , we tested our calculations 
 by artificially inverting the orientation of the diffusion tensors 
 of some of the structures and found very small changes in the relaxation 
 times as a result : a change of at most 0.25 ms in t1 and at most 0.02 ms in t2 .
these changes would not significantly impact the conclusions 
 of this manuscript and a detailed analysis is omitted here .
we 
 examined the five structures selected from the rdc - filter in light 
 of unbound tar rna structures generated using noe constraints ( but 
 not rdc constraints ) by one of the authors .
these 20 energy - minimized structures , termed the tar2013 series , 
 have  angles in the range 3859 , but the five 
 lowest energy structures cover a smaller range of 4553. 
 this range corresponds well to the bend angle of highest population 
 structure we have obtained ( i.e. ,  = 45 ) . in fact , the 
 lowest energy tar2013 structure has an  value of 202 
 and a  value of 225 , similar to the 
we previously published 
 two studies where the new methods developed to simulate relaxation 
 times were applied 
 to u38 relaxation data and used to select regions of interhelical 
 motional parameter space that fit the data .
the models consisted of 
 two - site jumps between the lowest energy structure of 1anr and structures 
 artificially modified from that structure to reflect changes in interhelical 
 orientation .
the closest approximation to a two - site jump in this 
 manuscript is found by considering only exchanges between pairs of 
 conformers within the three most populated structures ( the 45 
 structure , the 61 structure , and the 115 structure , which 
 are almost equal in population ) .
the exchange between the 45 
 and the 61 structures involves a bend angle modification of 
 16 , and a twist angle about the upper helix of 2. cross - checking 
 this parameter set against the results of applying the general rate 
 theory , we find that the u38 data was fit by a two - site 
 jump model with a twist of 0 and bend angle between 5 
 and 12 ( among other possible models ) .
thus , the 45 and 
 61 structures fit the profiles of two structures selected previously 
 on the basis of the u38 data alone .
exchanges between either of these 
 structures and the 115 structure , however , are of a magnitude 
 not simulated in the previous studies .
we can also make a few 
 basic comparisons to the results of dayie 
 et al .
, even though the authors consider 
 the hiv-2 tar rna , as both molecules consist of a helix  bulge  helix 
 motif . given our focus on the c6 and c8 atoms in the current work , 

we first observe a similarity in the fact that in the relaxation data 
 used in this work ( from bardaro et al . )
the helical residues seem to have similar c ( c6 and 
 c8 ) t1 relaxation times , a fact observed 
 in the work of dayie et al . as well
however , the u23 c t1 is nearly half the value 
 of the corresponding u25 time in dayie et al . , whereas bardaro et 
 al .
report values much closer in magnitude . also , the magnitudes of 
 the t1 and t1 relaxation times in the two papers are different .
we mention these 
 facts to bring up three relevant considerations in interpreting the 
 results of the two different sets of experiments : ( a ) the obvious 
 difference made by the presence of only two residues in the bulge 
 in hiv-2 rna versus three in the hiv-1 rna ; ( b ) the fact that the 
 relaxation times of dayie et al .
include a component from the c5c6 
 dipolar coupling for pyrimidines , whereas this coupling is explicitly 
 suppressed in bardaro et al .
( see shajani and varani ) ; ( c ) the use of a model - free analysis in dayie et al . 
 versus the se and gr methods used herein .
notwithstanding these caveats , 
 the common results we can extract are that both papers observe significant 
 flexibility in the u23 and u25 residues and rigid , slow motions in 
 the helices . to examine the extent to which our approach 
 matches results obtained 
 by the md approach
, we compare our results to those of salmon et al . , where the authors describe the selection of 
 conformers from 8.2 s md simulations of the tar rna using pf1 
 phage - aligned rdc data sets .
it was reported that the best fit to 
 rdc data is obtained with a set of 20 conformations selected from 
 the full md ensemble .
though it is not possible to compare the absolute 
 values of the interhelical bend and twist angles due to differing 
 methods of characterizing the helices and their relative orientations , 
 we can compare the spans of the angles reported for their ensemble 
 to those in ours .
ensemble of 
 20 span 88 ( from 3 to 91 ) whereas those in our ensemble 
 span 87 , the rotations of the upper helix about the lower helical 
 axis span 191 in their ensemble whereas those in our set span 
 215 , and the rotations of the upper helix about its own symmetry 
 axis span 224 in their ensemble and 210 in ours .
thus 
 the span of angles obtained by the two approaches are in excellent 
 agreement .
moreover , the full , prefilter ensembles in both papers 
 show correlations between the  and  twist angles .
a 
 more fine - grained comparison relates to the behavior of individual 
 residues .
we have already mentioned that the a22-u40 base pair is 
 often found to be open among the full set of 500 structures .
we also 
 find that , among our five rdc - filtered structures , four ( the 45 , 
 76 , 115 , and 132 structures ) lack a a22-u40 base 
 pair .
however , the g26-c39 base pair is maintained in all five of 
 these structures .
find a similar asymmetry between the 
 a22-u40 base pair and the g26-c39 base pairs in their rdc - selected 
 ensemble , with the former adopting a broader conformational distribution 
 and the latter being in an a - form helix - like conformation .
report the occurrence of three clusters within their rdc - selected 
 ensemble : a 66% population cluster with a22 stacked on u23 , a 19% 
 population cluster with u23 flipped out , and a 15% cluster with paired 
 u25-u40 and unpaired u23 and c24 . with regard to the third cluster , 
 nearly 30% of the 500 energy - minimized structures used in our analysis 
 show a u25-u40 pair , with three of the rdc - filtered structures ( the 
 45 , 76 , and 132 structures ) included in this list 
 as well .
the 115 structure simply lacks the a22-u40 pair and 
 does not have any alternative pairings of either residue .
stated that the u25-u40 pair is predicted to be the second most 
 energetically favorable bulge conformation in mc - fold .
a visual inspection 
 of our structures shows the following behavior for the bulge:(1)the 45 structure 
 has u23 flipped 
 into the interhelical space but not stacked , c24 is flipped in but 
 not stacked and u25 is paired with u40.(2)the 61 structure has u23 stacked 
 on a22 and has c24 and u25 flipped out.(3)the 76 structure has u23 flipped 
 out , c24 stacked with u25 , and u25 paired with u40.(4)the 115 structure has u23 flipped 
 out , c24 flipped in and stacked close to u25 , and u25 flipped in and 
 stacked only with c24.(5)the 132 structure has u23 and 
 c24 flipped in but not stacked and has u25 paired with u40 .
the 45 structure 
 has u23 flipped 
 into the interhelical space but not stacked , c24 is flipped in but 
 not stacked and u25 is paired with u40 .
the 61 structure has u23 stacked 
 on a22 and has c24 and u25 flipped out .
the 76 structure has u23 flipped 
 out , c24 stacked with u25 , and u25 paired with u40 .
the 115 structure has u23 flipped 
 out , c24 flipped in and stacked close to u25 , and u25 flipped in and 
 stacked only with c24 .
the 132 structure has u23 and 
 c24 flipped in but not stacked and has u25 paired with u40 .
thus , only one of the structures in our 
 ensemble ( the 61 
 structure with a 27% population ) has a significant a22-u23 stacking 
 interaction and two have u23 flipped out ( 39% total population ) , a 
 clear deviation from the results of salmon et al . , suggesting that 
 the conformational variability of this region is more than can be 
 captured by a small number of sampled structures .
a solution to this 
 problem is to generate energy - minimized structures where the interhelical 
 orientation is fixed ( or nearly so ) and the bulge flexibility is evaluated 
 under the constraint of fixed end points .
such an analysis would establish 
 the inherent conformational flexibility of the bulge . to cross - check the 
 results of the program hydronmr , we recalculated 
 the diffusion tensors using the program best , which tessellates the solvent - accessible surface of the molecule 
 and calculates the various diffusion properties using a finite element 
 analysis .
the eigenvalues 
 of the rotational diffusion tensors ( in ascending order ) of the two 
 programs are compared in table 4 .
we find that the diffusion 
 tensor eigenvalues as found by best 
 were uniformly smaller than those found by hydronmr , indicating that 
 hydronmr underestimated the hydration effect relative to best .
we subsequently 
 calculated the relaxation times using the slow exchange formalism 
 and found a t1 shift of at most 18 ms 
 and a t2 shift of at most 1.2 ms , corresponding 
 to a shift of about 4% of the experimental values for both relaxation 
 times .
this may modify the choice of parameters described above , but 
 we believe that the impact will not be substantial .
we have carried out a characterization of the 
 essential dynamics 
 of the tar rna molecule using techniques with time scale sensitivities 
 ranging from subnanosecond ( solid - state and solution relaxation times ) 
 to millisecond ( rdcs ) .
we have been able to capture the long - time 
 scale behavior of the conformational exchange processes that characterize 
 this molecule and fit experimental relaxation times very well , with 
 exchanges between discrete conformers occurring at time scales longer 
 than 1 s . the similarities of results of this method with those 
 of extended md simulations provide independent corroboration of our 
 conformational analysis .
further computational explorations and sample - size 
 increases will enhance the results obtained by this methodology .
 complex rna structures are constructed 
 from helical segments connected 
 by flexible loops that move spontaneously and in response to binding 
 of small molecule ligands and proteins . 
 understanding the conformational 
 variability of rna requires the characterization of the coupled time 
 evolution of interconnected flexible domains . to elucidate the collective 
 molecular motions and explore the conformational landscape of the 
 hiv-1 tar rna 
 , we describe a new methodology that utilizes energy - minimized 
 structures generated by the program  fragment assembly of rna 
 with full - atom refinement ( farfar )  . 
 we apply structural filters 
 in the form of experimental residual dipolar couplings ( rdcs ) to select 
 a subset of discrete energy - minimized conformers and carry out principal 
 component analyses ( pca ) to corroborate the choice of the filtered 
 subset . 
 we use this subset of structures to calculate solution t1 and t1 relaxation times for 13c spins in multiple residues in different domains of the molecule 
 using two simulation protocols that we previously published . 
 we match 
 the experimental t1 times to within 2% and the t1 times to within less than 10% for helical residues . 
 these results 
 introduce a protocol to construct viable dynamic trajectories for 
 rna molecules that accord well with experimental nmr data and support 
 the notion that the motions of the helical portions of this small 
 rna can be described by a relatively small number of discrete conformations 
 exchanging over time scales longer than 1 s . 
Introduction
Theoretical
and Computational Methods
Results
Discussion
Conclusions