Constructing RNA dynamical ensembles by combining MD and motionally decoupled NMR RDCs: new insights into RNA dynamics and adaptive ligand recognition.

We describe a strategy for constructing atomic resolution dynamical ensembles of RNA molecules, spanning up to millisecond timescales, that combines molecular dynamics (MD) simulations with NMR residual dipolar couplings (RDC) measured in elongated RNA. The ensembles are generated via a Monte Carlo procedure by selecting snap-shot from an MD trajectory that reproduce experimentally measured RDCs. Using this approach, we construct ensembles for two variants of the transactivation response element (TAR) containing three (HIV-1) and two (HIV-2) nucleotide bulges. The HIV-1 TAR ensemble reveals significant mobility in bulge residues C24 and U25 and to a lesser extent U23 and neighboring helical residue A22 that give rise to large amplitude spatially correlated twisting and bending helical motions. Omission of bulge residue C24 in HIV-2 TAR leads to a significant reduction in both the local mobility in and around the bulge and amplitude of inter-helical bending motions. In contrast, twisting motions of the helices remain comparable in amplitude to HIV-1 TAR and spatial correlations between them increase significantly. Comparison of the HIV-1 TAR dynamical ensemble and ligand bound TAR conformations reveals that several features of the binding pocket and global conformation are dynamically preformed, providing support for adaptive recognition via a 'conformational selection' type mechanism.

INTRODUCTION

Many non-coding RNA molecules (ncRNAs) perform their biological functions by undergoing large conformational changes in response to specific cellular signals including the recognition of proteins, nucleic acids, metal ions, metabolites, vitamins, changes in temperature and even RNA biosynthesis itself (1–5). These conformational transitions guide RNA folding during co-transcriptional folding; provide the molecular basis for sensing and signaling transactions that allow riboswitches to regulate gene expression in response to changes in environmental conditions; allow ribozymes to dynamically meet the diverse structural requirements associated with their multi-step catalytic cycles; and enable complex ribonucleoproteins to assemble in a hierarchical and sequentially ordered manner.

Although it is clear that many ncRNAs undergo large changes in structure in order to carry out their function, the mechanism by which these conformational transitions occur remains poorly understood. A central question is whether cellular factors such as proteins and ligands act catalytically to induce the RNA conformational change via ‘induced fit’, or they select and bind distinct RNA conformers from a pre-existing dynamical ensemble via ‘conformational selection’ (6–8). Insights into such mechanistic questions have been impeded by lack of biophysical techniques that allow the 3D visualization of intrinsic RNA dynamics over biologically relevant timescales. The atomic resolution characterization of dynamics in complex biomolecules is currently a major challenge in structural biology and biophysics. NMR spectroscopy is one of the most powerful techniques for characterizing dynamics uniquely providing comprehensive information regarding the amplitude, timescale and—in favorable cases—direction of motions with site-specific resolution (9–11). However, even with abundant measurements that can be made with the use of NMR, the total number of observables still pale in comparison to the total number of parameters needed to fully describe dynamics. Molecular dynamic (MD) simulations provide an all-atom description of dynamics; however, force fields remain to be thoroughly validated particularly for nucleic acids and simulation timescales remain limited to ∼100 ns (12,13).

Because they are complementary on the spatial and temporal scales, the limitations inherent to NMR and MD could in principle be overcome by combining the two techniques; MD can fill the shortage in NMR data and NMR can provide a means for validating and potentially correcting force fields and accelerate MD conformational sampling to millisecond timescales. Several studies have emerged in which MD and NMR are used in concert in studies of nucleic acid dynamics (5,14–17). The two techniques can also be integrated to yield a unified view of structural dynamics. While such combined NMR/MD approaches have successfully been used in studies of protein dynamics (18–21), extension to RNA can prove very difficult. This is because unlike globular proteins, overall motions in highly flexible RNAs can be strongly coupled to internal motions making it difficult, if not impossible, to predict NMR data from an MD trajectory (22–24).

Recently, a domain-elongation strategy was introduced for decoupling internal and overall motions in RNA (25,26). In this approach, a target helix is elongated so that the overall motion is slowed down relative to internal motions and rendered less sensitive to internal fluctuations in other parts of the molecule. The elongation also simplifies analysis of NMR spin relaxation (27) and residual dipolar coupling (RDC) (28,29) data because it predefines the overall diffusion or alignment tensor to be axially symmetric with principal axis oriented nearly parallel to the elongated helix axis. As we show in what follows, this makes elongated RNAs ideally suited for computing NMR observables from a given MD trajectory of a corresponding non-elongated RNA in which snapshots are aligned by superimposing the reference elongated helix.

Here, we report implementation of a general strategy that combines MD simulations and NMR RDCs measured in elongated RNAs for constructing atomic resolution dynamical ensembles with timescale sensitivity extending up to milliseconds. The dynamical ensemble of the HIV-1 and HIV-2 transactivation response element (TAR) (Figure 1a ) constructed in this manner provides new insights into the bulge-length dependence of RNA dynamics and the mechanism of conformational adaptation upon target recognition.

Figure 1.

SAS analysis of HIV-1 E-TAR RDCs. (a) Secondary structure of HIV-1 TAR with helix I highlighted in red, helix II in green and tri-nucleotide bulge in orange. HIV-2 TAR lacks bulge residue C24. (b and c) Plots of experimental RDCs versus values computed from the 80 ns MD trajectory for (b) EI-TAR and (c) EII-TAR. Data for helix I, helix II, and bulge, are shown in red, green and orange, respectively. Also shown is the root-mean-square-deviation (rmsd) and correlation coefficient (R). (d) RMSD (Hz) between calculated and experimental RDCs as a function of N, following a SAS analysis using both EI and EII-TAR HIV-1 RDCs. (e–g) Plots of experimental RDCs versus values calculated from the N = 20 SAS ensemble using (e) EI-TAR, (f) EII-TAR and (g) EI-TAR+EII-TAR RDCs.

MATERIALS AND METHODS

Molecular dynamics simulations

Simulations of wild-type HIV-1 and HIV-2 TAR were performed using the CHARMM package (30) with force field parameter set 27 (31). Model 3 of the unbound NMR structure of HIV-1 TAR (PDB 1ANR) (32) was used as starting coordinates for simulations of HIV-1 TAR. For HIV-2 TAR, starting coordinates were obtained by removing argininamide from model 1 of the argininamide-TAR NMR structure (PDB 1AKX) (33). The RNA was charge-neutralized using sodium counter ions and solvated in a 33–35 Å sphere of TIP3 water (34). A spherical boundary potential was applied to maintain the density of water around the RNA site (35). The 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/Å for 100 ps, after which constraints were removed and the system equilibrated for 1 ns. A Nosé–Hoover thermostat (36,37) was used to maintain a constant temperature of 300 K throughout the simulation, with a 1 fs time-step and a coupling constant of 50 ps−1. Fifty distinct trajectories were initiated from this equilibrated structure, by assigning different initial velocities. It has previously been demonstrated that this technique can be used to enhance conformational sampling relative to a single trajectory of the same total duration (38,39). The first 0.5 ns of each trajectory was discarded and next 1.6 ns used for analysis. Conformations from each of the 50, 1.6 ns, trajectories were pooled to give a total effective simulation time of 80 ns. These 80 000 conformations were used as a structural pool for the ‘selection’ phase of the select-and-sample strategy (SAS, see below). The same protocol used for HIV-1 TAR, was used to prepare, equilibrate and produce the MD trajectories for HIV-2 TAR, except that a 33 Å sphere of TIP3 water was used to solvate the system.

RDC-based SAS

In the original implementation, Chen and coworkers used S2 NH spin relaxation order parameters to select structural ensembles for a series of protein systems (19). Here, we adapt this SAS approach to allow use of RDC data which probe the orientational dynamics of individual bond vectors over timescales extending up to milliseconds. In the SAS approach, one generates a set of conformations for the system of interest, in which we use MD to generate such structures. One then searches for an N-membered subset of structures that minimizes a cost function. Specifically, an N-membered subset of structures is randomly selected from a total pool of M structures and an initial χ2-value is evaluated using Equation 1. Next, one of the N-membered structures is randomly chosen and replaced, by a random structure from the remaining M–N conformational pool. The ‘move’ from step k to k + 1 is then accepted if χ2(1) < χ2(). If χ2(+1) > χ2(), the move is accepted with a probability P = exp((χ( −χ(+1))/Teff), where Teff is an effective temperature that is linearly decreased in a simulated annealing scheme. The cost function is where and are the calculated and measured RDCs, respectively, L is the total number bond vectors, and is calculated using, where γi is the gyromagnetic ratio of the i-th nucleus, rij is the bond length, θ is the angle between the inter-nuclear bond vector and the external magnetic field and the angular brackets denote a time-average over all sampled orientations. The time-averaged angular term can be expressed in terms of the time-independent orientation of an internuclear vector relative to an arbitrary frame and the five order tensor elements (S), (40,41). where α is the angle between the internuclear vector and the n-th axis of the arbitrary frame. Equation 3 assumes that internal motions do not affect overall alignment of the molecule. This assumption can break down in RNA because collective motions of helical domains can lead to large changes in the overall shape and thus overall alignment of the molecule (22,23,42). This can make it impossible to accurately compute Dcal for a given ensemble and thus to use RDCs in selecting conformers form an MD trajectory. Domain-elongation allows one to overcome this problem by rendering the overall shape of the molecule far less sensitive to internal motions (25,26). Elongation also predefines the overall order tensor to be axially determined with principal direction oriented approximately parallel to axis of the elongated helix (26). This overall order tensor can be conveniently determined experimentally using RDCs measured in the elongated helix (26).

TAR ensembles by RDC-based SAS

The SAS RDC approach was implemented using in-house programs written in C++ available from the authors upon request. The previously reported (26) 47 and 35 one-bond base and sugar C-H RDCs measured in EI-TAR and EII-TAR, respectively were used independently or in combination to select an HIV-1 TAR structural ensemble from a pool of 80 000 conformers derived from the combined 80 ns MD trajectory. For HIV-2 TAR, 35 RDCs measured in EI-TAR were used. Note that although RDCs were measured on a TAR construct in which the apical loop was replaced by a UUCG loop, a detailed NMR study recently showed that this apical loop replacement does not affect inter-helical motions or local motions at the TAR bulge (43). When used in combination, both the EI-TAR and EII-TAR sets of RDCs were used (L = 47 + 35 = 82) in computing the χ2 in Equation 1. In all cases, the overall alignment tensor was assumed to be axially symmetric (η = 0) with principal direction (Szz) oriented parallel to the elongated helix axis as computed using the program CURVES (44,45). The magnitude of Szz was obtained from an order tensor analysis as described previously (26).

Each RDC–SAS run was initiated from N randomly selected conformers. A Monte Carlo (MC) simulated annealing scheme was then used to minimize the cost in Equation 1 as described above. Simulations were started at a high-effective temperature, where the MC acceptance probability was high (∼0.99), and slowly decreased until the MC acceptance probability was ∼10−5. At a given effective temperature 105 MC steps were carried out. The effective temperature was then decreased, with T+1 = 0.9T. The same protocol was used for HIV-2 TAR, except that only 35 RDCs measured in the helix I elongated sample were used in the SAS analysis. Here, there was less motivation to acquire an additional set of helix II elongated RDCs given the near linear alignment of the two helices which renders the two sets of data degenerate.

Analysis of the SAS-derived TAR ensemble

The base angles, buckle (κ), opening (σ), propeller twist (ω) and twist (Ω), were calculated for eight non-terminal base pairs for each of the SAS selected TAR conformers using the program 3DNA (46). The inter-helical Euler angles αh, βh and γh were also computed for each member of the ensemble, using the lower helix as a reference, as previously described (47). Note that degenerate sets of Euler angles reflect the same inter-helical orientation (47) and care had to be taken in selecting among degenerate angles to avoid non-realistic correlations.

RESULTS

SAS analysis of HIV-1 E-TAR RDCs

In Figure 1, we compare the RDCs measured previously (26) in EI-TAR and EII-TAR with those predicted based on the combined 50 × 1.6 = 80 ns MD trajectory. Though some correlation is observed between the measured and predicted RDCs for both EI-TAR (Figure 1b) and EII-TAR (Figure 1c), the deviations (RMSD ranging between 13 and 16 Hz) remains substantially larger than the estimated RDC measurement uncertainty (∼4 Hz). The MD trajectory does not systematically under- or overestimate the RDCs measured throughout the RNA, indicating that it does not significantly over- or underestimate the amplitude of motions present. This is the case even though the RDC timescale sensitivity to motions (less than millisecond) is greater than that of MD (∼80 ns). Thus, it appears that the motions in TAR saturate at the nanosecond timescales consistent with relaxation dispersion studies that provide no evidence for µs-ms motions in the stem-bulge-stem element of TAR (43).

We examined if SAS could be used to pull out sub-ensembles from the MD trajectory that satisfy the measured RDCs. We first conducted a series of SAS runs with various ensemble sizes of N = 5, 10, 20, 30, 40, 50 and 100. The root-mean-square-derivation (RMSD) between measured and predicted RDCs when combining the EI-TAR and EII-TAR RDCs is shown in Figure 1d as a function of N. Increasing the ensemble size beyond N = 20 did not lead to significant improvements in the fit and in fact a deterioration was observed for N > 20. This is likely due to sampling problems during the MC simulating annealing minimization of the cost function as the number of possible combinations increases steeply with N and the cost function exhibits, due to the frustration stemming from the underdetermined nature of the problem, a vast number local minima. An ensemble size of N = 20 was used in all subsequent SAS runs.

The RDCs calculated using a 20-member ensemble selected using the SAS approach exhibit a markedly improved fit to the EI-TAR (Figure 1e), EII-TAR (Figure 1f) and EI-TAR+EII-TAR (Figure 1g) RDCs as compared to those calculated from the entire MD trajectory. Importantly, sub-ensembles can be determined that simultaneously reproduce the EI-TAR and EII-TAR RDCs with an RMSD (4.8 Hz) that is comparable to the estimated experimental RDC uncertainty (∼4 Hz) (Figure 1g). The SAS selected conformers also did not lead to any steric collisions with the elongated helices for both EI-TAR and EII-TAR despite the fact that the helices were not actually elongated in the MD simulations (data not shown).

To examine the uniqueness of the selected TAR ensemble, we carried out 100 independent N = 20 SAS runs. Here, a maximum of 2000 unique snapshots can be selected. However, repeated runs resulted in selection of a narrow set of similar conformations. The total number of unique snapshots selected by the RDC-SAS optimization algorithm was 422, 48 and 121 for EI-TAR, EII-TAR and EI-TAR+EII-TAR RDCs, respectively. Thus, the RDC data favors selection of specific conformations from the available pool. As expected, the conformers selected did vary when changing the value N. However, the overall distribution of conformations remained similar as shown for example for the inter-helical orientation in Supplementary Figure S1.

Local motions in the helices and bulge

In Figure 2a , we plot (in black) the average and standard deviation for various base angles computed for the 121 unique TAR conformers obtained from 100 N = 20 SAS runs. For comparison, the mean value and SD for idealized A-form helical geometry computed from a statistical comparison of high-resolution X-ray structures is shown in red (48). For the majority of the residues, very good agreement is observed between the SAS ensemble and canonical values. Large deviations are however observed for the junctional A22-U40 base-pair, which in the MD simulation frequently deviates from a hydrogen bonded alignment. This is in excellent agreement with previous NMR data showing that while the junctional G26-C39 base-pair forms a detectable hydrogen bond, the A22-U40 base-pair is flexible and does not form the expected base-pair in HIV-1 TAR (49–51).

Figure 2.

Local motions in the HIV-1 TAR dynamical ensemble. (a) Shown in black are the mean values for the base opening (σ), buckle (κ), propeller twist (ω) and twist (Ω) angles and their SD calculated over 121 HIV-1 TAR conformers obtained from multiple N = 20 SAS runs. For comparison, shown in red are corresponding values for an idealized A-form helix as obtained from a statistical survey of high resolution X-ray structures (48). (b) Average conformation of the HIV-1 TAR bulge and neighboring base-pairs calculated from the 121-membered SAS ensemble. The bases of the bulge and flanking base pairs are color-coded based on the root mean square fluctuations (r.m.s.f.) calculated over the ensemble.

The bulge residues exhibit different levels of motions. As shown in Figure 2b, the highly conserved U23 bulge stacks onto A22 in the majority of the conformations, consistent with observation of NOE connectivity between A22 and U23. The two residues undergo limited motions consistent with previous 13C relaxation studies of dynamics in elongated HIV-1 TAR (52). Interestingly, select conformations exist in which U23 adopts a looped out conformation as observed in several ligand bound TAR structures (PDB ID# 1QD3, 1UTS, and 397D). In contrast, residues C24 and U25 predominantly exist in a looped out conformation and are significantly more flexible, again in agreement with the previous carbon relaxation studies (52).

Collective motions of helices

The relative orientation of two helices i and j can be defined using three inter-helical Euler angles that describe the twist angle αh about helix i, the inter-helix bend angle βh, and the twist angle γh about helix j (26,47). We computed these three inter-helical Euler angles for the 121 SAS selected TAR conformers and compared them to angles obtained for the entire MD trajectory. As shown in Figure 3a, the MD trajectory spans a large range of inter-helical angles. The SAS conformers are widely distributed across the MD trajectory (Figure 3a). For both the MD trajectory and SAS ensembles, significant correlation is observed between the three inter-helical angles, particularly between the twist angles αh and γh. Such spatial correlations were reported previously based on an a three-state rigid-body refinement of the TAR inter-helical orientation (26). As shown in Figure 3b, the conformations obtained by SAS sample orientations that are in very good agreement with those obtained by a previous three-state analysis of the RDCs (26). The only significant deviations are observed in the αh-γh plane for the near coaxial conformer (βh∼ −21°). This is not surprising given that for this near coaxial conformer, the RDCs measured in both EI-TAR and EII-TAR are highly insensitive to the twist angles αh and γh.

Figure 3.

Global inter-helical dynamics in the HIV-1 TAR dynamical ensemble. (a) Shown in gray are the inter-helical twist (αh and γh) and bend (βh) angles for 80 000 HIV-1 TAR conformers derived from an 80 ns MD trajectory (conformer selected every 1 ps). In blue are the corresponding SAS selected conformers. The correlation coefficient (R) is shown on individual planes. (b) Comparison of the SAS selected inter-helical angles and those derived previously (26) based on a three-state rigid body ensemble analysis of E-TAR RDCs.

We also examined if there exists any correlation between the geometry of base-pairs at the junction and the inter-helical bend angles. The only significant observation was anti-correlation (R ∼ −0.70) between the base-pair step angle (Ω) at the G26-C39 junctional base-pair and αh the twist angle about helix II. Thus, the unraveling of the G26-C39 base-pair may give rise to twisting motions around the axis of helix II.

Comparison with HIV-2 TAR motions

We used the SAS approach to analyze RDCs previously measured in HIV-2 EI-TAR (26) in which the bulge residue C24 is omitted. Poor agreement (RMSD = 15.1 Hz) was again observed between the measured EI-TAR RDCs and values computed using entire MD trajectory (Figure 4a ). By using SAS, we were able to find an N = 20 sub-ensemble that yields an RMSD of 1.7 Hz (Figure 4b and Supplementary Figure S2). Compared to HIV-1 TAR, repeated SAS runs resulted in selection of a larger number (276) of unique conformers for HIV-2 TAR. This could be attributed to a smaller RDC sensitivity to twisting motions both because the HIV-2 TAR structure is more linear and because only RDCs measured in EI-TAR was available for analysis.

Figure 4.

HIV-2 TAR dynamical ensemble. Plot of experimental HIV-2 EI-TAR RDCs versus values calculated from (a) an 80 ns MD trajectory of HIV-2 TAR and (b) following SAS selection of N = 20 conformers from the MD trajectory based on HIV2 EI-TAR RDCs. Coloring scheme is same as in Figure 1. (c) Shown in black are the mean base-pair angles and their SD calculated over a 279-membered HIV-2 TAR SAS ensemble. For comparison shown in green are corresponding angles for an idealized A-form helix as obtained from a statistical survey of high-resolution X-ray structures (48). (d) Average conformation of the bulge of HIV-2 TAR calculated from the 279-membered HIV-2 TAR SAS ensemble. Bases of the bulge and flanking base pairs are color-coded based on the root mean square fluctuations (r.m.s.f.) calculated over the ensemble. (e) Comparison of SAS selected HIV-1 (blue) and HIV-2 (red) inter-helical conformations. The correlation coefficient (R) is shown on individual planes.

Comparison of the SAS selected HIV-2 TAR conformers with those obtained for HIV-1 TAR (see Supplementary movies) revealed that reducing the length of the bulge in HIV-2 TAR led to a marked decrease in the local motions in the junctional A22-U40 base-pair (Figure 4c). In HIV-1, the SDs are approximately 30°, 45°, 13°, 16° for the base angles κ, σ, ω and Ω, respectively. In HIV-2, they reduce to 10°, 6°, 8° and 6°, respectively. In contrast, we observe significant static deviations in the opening angle for the G26-C39 junctional base pair. Likewise, a significant reduction is observed in the local dynamics of bulge residues U23 and C25 (Figure 4d). U23 is less flexible and forms more stable stacking interactions on an also less flexible A22. The root mean square fluctuations (r.m.s.f.) of the atomic positions of U23 and U25 decrease from 1.89 Å and 3.49 Å, respectively in HIV-1 TAR to 1.45 Å and 1.49 Å, respectively in HIV-2 TAR.

The reduction in the local motions in and around the bulge linker is, as expected, accompanied by a reduction in the inter-helical motional amplitudes, as shown in Figure 4e. Such a reduction is clearly observed for the inter-helical bending which decreases in standard deviation from ∼33° to ∼12°. This is in agreement with an order tensor analysis of RDCs, which reported a reduction in the ϑint (which ranges between 0 and 1 for maximum and minimum inter-helical motional amplitudes, respectively) from 0.45 ± 0.05 to 0.77 ± 0.04 (26). As expected, the HIV-2 conformers cluster more tightly around more linear (βh ∼0°) conformations. While we do not observe a significant reduction in the amplitude of twisting motions about the two helices (αh and γh) in HIV-2, the RDC sensitivity to these angles is diminished in HIV-2 TAR both because the structure is on average more linear and because RDCs were only measured on the domain I elongated construct. Nevertheless, we observe a stronger correlation between the twisting motions indicating that the αh and γh correlations originates in part from the steric drag one helix exerts on the other.

Dynamics and adaptive recognition of ligands

Numerous studies have shown that HIV-1 TAR undergoes large conformational rearrangements that allow binding of diverse targets in and around the bulge, including peptide derivatives of its cognate protein Tat (53–55), divalent ions (56), and five chemically distinct small molecules (57–60). The three inter-helical conformers obtained by a three-state ensemble analysis of RDCs measured in unbound TAR revealed a global inter-helical motional trajectory that encapsulated many of the ligand bound conformations (26). As shown in Figure 5a , the SAS selected inter-helical conformers trace orientations similar to the ligand-bound TAR conformations, supporting the notion that unbound TAR can dynamically access its ligand-bound global conformations. This can also be seen in Figure 5b in which we compare the best matching ligand-bound TAR conformers and SAS-selected conformers as determined by superimposing all heavy atoms excluding the flexible terminal base-pairs (G17-C45) and the apical loop.

Figure 5.

Comparison of the SAS derived HIV-1 TAR dynamical ensemble and ligand bound TAR conformations. (a) Comparison of the global inter-helical angles. Shown in blue are the SAS selected angles and in gray seven distinct ligand bound TAR structures. (b and c) Comparison of the b, global and c, local structure of SAS TAR conformers and seven distinct ligand bound TAR conformations (1QD3, 1UUI, 1UTS, 1UUD, 1ARJ, 1LVJ, and 397D). Shown are the pairs yielding the lowest RMSD fit when superimposing b, all heavy atoms excluding terminal base-pairs G17-C45 and the apical loop and c, all heavy atoms in the bulge and immediately adjacent base-pairs. Every model in the ligand bound NMR ensembles was used in the superposition. The corresponding ligand is colored yellow.

The SAS dynamical ensembles also allowed us to examine to what extent are local features of the ligand binding pocket dynamically preformed in the absence of ligands. In Figure 5c, we compare the best matching ligand-bound TAR conformers and SAS-selected conformers as determined by superimposing heavy atoms in the bulge and immediately neighboring base-pairs. In general, the largest deviations are seen for the highly flexible residues C24 and U25, which are also known to be flexible in many of the ligand bound TAR conformations (26,61,62). Overall, these results suggest that local dynamics in an around the TAR bulge likely facilitate formation of the ligand binding pocket.

DISCUSSION

RNA elongation provides a reliable approach for computing time-averaged RDCs in highly flexible RNAs providing a basis for integration with MD as we showed here using the SAS approach. That the measured RDCs in both HIV-1 and HIV-2 can be satisfied using conformers selected from the MD trajectory suggests that neither HIV-1 or HIV-2 TAR undergo significant μs-ms motions that are not sampled in the MD trajectory. This is consistent with relaxation dispersion NMR studies of HIV-1 TAR which provide no evidence for μs-ms motions in and around the bulge (43). In further support of this idea comes the observation that the first few low-frequency eigenmodes of quasiharmonic motion obtained by principal component analysis of the 80 000 snapshots capture the dynamics as described by the SAS ensemble (Supplementary Figure S3).

However, it should also be noted that the SAS conformers represent a ‘discrete’ approximation to what is more likely a continuous complex distribution of many more conformations. The SAS conformers can be thought of as discrete points along the configuration space that may help define salient features of the motional trajectory (63). It should be kept in mind that averaging over such a discrete number of conformers leads to efficient averaging of the RDCs. A correspondingly larger continuous distribution of many more conformations will likely be required to accomplish the same level of motional averaging. It is therefore not surprising that the MD trajectory evolves outside the envelope defined by the SAS conformers and that the SAS envelope increases with the size of the ensemble (Supplementary Figure S1). Thus, the most likely source of discrepancy between the measured RDCs and the MD simulation is the assignment of the relative weights to the various conformers. RDC studies of partially unfolded proteins have emphasized the exquisite sensitivity of RDCs to the underlying conformational distribution (64).

The comparison of the dynamical ensembles generated for HIV-1 and HIV-2 TAR provided fundamental new insights into the dependence of RNA dynamics on the bulge length. As would be expected, reducing the length of the HIV-1 trinucleotide bulge by a single nucleotide (HIV-2) led to a significant reduction in the local motions in and around the bulge, as well as global inter-helical motions (see Supplementary movies). In both HIV-1 and HIV-2, we observe spatial correlations between twisting motions about individual helices, and to a lesser extent between twisting and bending as first reported based on the three-state ensemble analysis of the RDCs (26). Importantly, the spatial correlations between the twisting motions increase significantly in HIV-2 TAR (Figure 4e). This is expected if one were to consider limiting cases for the bulge length; at the limit of not having a bulge linker, the twisting dynamics becomes perfectly correlated, whereas for an infinitely long bulge, one would expect little correlation.

The HIV-1 TAR dynamical ensemble allowed us to directly examine if unbound TAR dynamically samples diverse ligand-bound conformations that have been reported to date. Our results suggest that many of the key features of the ligand bound TAR conformations, including the global inter-helical orientation and local aspects of the ligand-binding pocket appear to be dynamically preformed in unbound TAR. However, a detailed comparison of the SAS conformations and the ligand bound TAR structures remains complicated by a number of factors. As mentioned above, the SAS conformers only represent an approximate discrete state ensemble to what is likely a more continuous conformational distribution. One also has to consider the uncertainty in ligand bound TAR structures, both due to experimental imprecision and because the ligand bound states may be flexible themselves. Previous NMR studies reveal significant mobility in the TAR bulge, especially residues C24 and U25 when in complex with argininamide (49,65) and Mg2+ (62,66). The ACP bound TAR structure also exhibits a large degree of inter-helical motions (61). Notwithstanding these complications, our results suggest that local dynamics in and around the bulge together with global motions of helices can drive much of the conformational adaptation required to bind to different ligand targets.

Our results also suggest that some of the TAR conformational changes require ligand binding to occur efficiently. For example, we do not observe the key U38-A27•U23 base-triple in either of the HIV-1 or HIV-2 SAS conformers that is known to form in the TAR-argininamide complex (53). In this regard, it is interesting to note that larger deviations between the SAS and ligand bound conformations are generally observed for the weaker binding ligands and vice versa. It is possible that the weaker binding ligands expend a greater fraction of the binding energy changing the TAR conformation. Additional studies are needed to shed light on these key energetic questions.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online.

FUNDING

NSF CAREER Award (CHE-0548047 to I.A.); the National Institutes of Health (RO1 AI066975-01 to H.M.A.); National Science Foundation Graduate Research Fellowship Program (to A.C.S. and A.T.F.).

Conflict of interest statement. None declared.