The incorporation of the experimental SAXS data I exp ( q ) to MD simulations (SAXS-driven MD) was achieved by adding an additional energy term to the Hamiltonian E _hybrid = E _FF + E _SAXS , where E _FF is the energy from the MD force field, whereas E _SAXS represents the energetic penalty if the real-time computed amplitude I c o m p ( q ) deviates from the experimentally measured scattering curve I exp ( q ) . The real-time SAXS curves are calculated on the fly from the simulation using the method described in by Chen and Hub (2015) and He et al. (2021). Typically, the energy penalty between the simulated conformations and the experimental data is evaluated by an uncertainty-weighted coupling term, E S A X S R , I exp , t = α t k c k B T n q ∑ i = 1 n q I c o m q i , R , t − I exp q i 2 σ 2 q i , (4) where σ ² (q _i ) is the experimental error given as σ 2 ( q i ) = σ exp 2 ( q i ) + σ c o m 2 ( q i ) + σ b u f f e r 2 ( q i ) that estimates the total error as the sum of experimental error, statistical error from computed curves, and systematic error from the uncertainty of the buffer density, respectively. I c o m ( q i , R , t ) is the real-time SAXS curve back-calculated from the simulation, which is the temporal-weighted average within a memory time of t, whereas n _q is the number of intensity points spanning the specific range of scattering vector q. The coefficient k _c is a constant that adjusts the weight of the SAXS potential E _SAXS , as compared to the force field E _FF term. k _B and T denote the Boltzmann constant and temperature, respectively. Lastly, the parameter α(t) is a time-dependent function that allows a gradual introduction of the coupling potential at the start of the simulation.

The coupling protocol in the previous section assumed that the ensemble could be dominated by a single, well-defined structure. This is often a useful approximation when the RNA molecule of interest has a known and well-defined structure such as double helix topologies of RNA and DNAs (He et al., 2021). However, some RNA molecules show structural heterogeneity. The resulted experimental signal, in our case, the SAXS profile, comes from a heterogeneous ensemble of distinct conformations. As a result, a conformational pool around one unique structure may be impossible to successfully satisfy the experimental restraints. The maximum entropy principle provides a unique and a minimally biased choice of restraint potential that can be used to bias the system for heterogeneous ensemble such that the experimental restraints are satisfied on average by introducing minimal unphysical artifacts (Pitera and Chodera, 2012). Introducing an experiment-derived biasing energy (E _SAXS ) into the simulation gives a biased ensemble p _b (R) that represents the data that, according to the maximum entropy principle, are assumed to be as similar as possible to the original, unbiased ensemble p _ub (R). The similarity could be mathematically evaluated by Kullback–Leibler divergence D K L p b | p u b = ∫ p b R ln p b R p u b R d R . (5)

Flexible structures are expected to adopt heterogeneous conformations; hence, the biased ensemble distribution p _b (R) obtained by restraining a single simulation to the data (Eq. 4) deviates from the original ensemble distribution p _ub (R), resulting in a high D _KL , and that goes against the maximum entropy principle. To minimize the D _KL , additional biasing terms of the back-calculated experimental data and the weight of each conformation are determined by Lagrange multipliers whose value is computed to enforce agreement between experiment and simulation.

To avoid the complications of determining the Lagrange multipliers, the replica-averaged modeling has been proposed (Best and Vendruscolo, 2004; Lindorff–Larsen et al., 2005). In this approach, a set of replicas of the system is simulated in parallel. The parallel-replica simulations with harmonic constraint to the ensemble-averaged data show to comply with the maximum entropy principle (parallel-replica SAXS-driven MD) (Pitera and Chodera, 2012; Roux and Weare, 2013). In this scheme, the simulated scattering curves are back-calculated by averaging the parallel replicas and giving more weight to the recent conformations (Chen and Hub, 2014) I ̄ c o m q i , R 1 , … , R N , t = N − 1 ∑ β = 1 N I c o m q i , R β , t , (6) where the replica index is denoted as β, so Eq. 4 of experiment-derived energy penalty is then modified accordingly E S A X S R 1 , … , R N , I exp , t = α t k c N Ω k B T n q ∑ i = 1 n q I ̄ c o m q i , R 1 , … , R N , t − I exp q i 2 σ 2 q i (7)

Following the works of Hummer and Köfinger (2015) and Hermann and Hub (2019), setting the exponent Ω in Eq. 7 to Ω = 1 ensures that the optimal ensemble is recovered as N → ∞. Finally, it is to be noted that Eq. 4 and Eq. 7 share key parameters.

We study the SAM-I riboswitch aptamer at three solvent conditions. That is, (i) the SAM-bound and (ii) the SAM-free RNA molecule in the presence of Mg²⁺ ions (denoted henceforth as s1 and s2, respectively), and (iii) the Mg²⁺- and SAM-free state of the riboswitch (s3). We used for all these conditions the crystal structure of the SAM-bound complex (PDB entry: 2GIS) as the initial structure due to its good agreement with the reported SAXS data (Stoddard et al., 2010). The ligand was removed for the study of s2 and s3. Each of these starting coordinates was placed into a triclinic box with sufficient space to ensure structural changes of the RNA molecule during MD simulation toward the SAXS-based structural information. In detail, we selected a distance between the solute and the box of 1.2 and 4 nm for the conditions s1 and s2-3, respectively. The SAM-I riboswitch aptamer was solvated with water and ions to mimic the experimental condition of 7.6 mM MgCl₂ and 150 mM NaCl. Additional 93 Na⁺ ions were added to satisfy charge neutralization. All these mono- and divalent ions were added to the simulation box by replacing water molecules. Energy minimization was employed to refine the structures using GROMACS 5.0.5. We used the χOL3 force field (Zgarbová et al., 2011) to represent the RNA in combination with the TIP3P water model (Jorgensen et al., 1983). The force field was extended to include the topological information for SAM using the previous work (Hayes et al., 2012). The experimental SAXS curves show that the structure of the apo-state dramatically changes in the absence of Mg ²⁺ ions (s3). Hence, we used different starting structures. These starting structures were prepared by partially unfolding the RNA by separating the helices either through heating or pulling. For the heating simulations, the crystal structure was equilibrated in a 150 mM aqueous solution by running 5-ns-long simulation in an isothermal–isobaric (NPT) ensemble at 300 K, followed by a 20-ns-long canonical (NVT) ensemble at 350 K. To disrupt the specific tertiary contacts leading to the compact structures, we also applied force. After a 5-ns NPT simulation run at 300 K, the helices were pulled apart during a 20-ns NVT simulation with a pulling rate of 0.001 nm/ps and a force constant of 1,000 kJ/(mol × nm²). The pulling direction was defined as the x-component of the center of mass distance between P1 (nucleotides A6 and G89) and P3 (nucleotides U57 and C47). The two 20-ns-long trajectories were analyzed using the GROMOS algorithm with a cutoff of 0.2. The SAXS profile was calculated for each resulting cluster. The SAXS curve of seven structures was found to be close to the reported experimental SAXS curve for the Mg²⁺-free apo-state. Each of this structure was then energy minimized analogous to the s2 structure.

The starting structure of the RNA HJH construct, including two 12 base-paired RNA helices linked by a short chain of rU₅, was built using Nucleic Acid Builder (NAB) (Case et al., 2005) by assuming an A-form geometry for the helix. The long backbone has a mixed sequence as CCC UAU ACU CCC UUU UUC CUC CUA AUC GC. The construct was placed in a cubic box by aligning the long axis parallel to the z-axis and energy minimized using GROMACS 5.0.5 with the same force parameters as the aforementioned ligand-free riboswitch system. The simulation box along the z axis was extended by about 20 Å at both ends of the RNA HJH to ensure that the molecule does not interact with its image. The resulting cubic box had a side length of 13.0 nm. The system was then solvated with TIP3P water and ions were added to match the experimental conditions. We studied the HJH system in the absence of K⁺ ions and to be consistent with experiment, in the presence of 50 mM, 100 mM, or 500 mM K⁺ concentration.

All MD simulations of the riboswitch were carried out using the GROMACS 5.0.5 suite (Hess et al., 2008) with the imposed periodic boundary conditions in all directions. Possible bad contacts resulting from the random water and ion placement were removed by a 50,000-step long steep descent minimization. The minimized structures were then equilibrated for volume and solvent as follows: 5-ns-long simulation in isothermal–isobaric ensemble was employed with a constant temperature of 300 K and pressure of 1 bar ensured by using Berendsen thermostat (Berendsen et al., 1984) and Parrinello–Rahman barostat (Parrinello and Rahman, 1981), respectively. During this step, the heavy atoms of the RNA were restrained to their initial positions during the simulation, whereas ions and water were allowed to move freely. For this, harmonic restraints were applied to the RNA atoms with a force constant of 1,000 kJ/nm². Particle mesh Ewald (PME) (Darden et al., 1993) summation was used to compute long-range electrostatic interaction with a grid spacing of 0.16 nm and an interpolation of order 4. The real space distance cutoff (for electrostatics and van der Waals energies) was set to 10-Å covalent bonds of the water and nucleic acid, which were constrained to their equilibrium geometries using SETTLE (Miyamoto and Kollman, 1992) and LINCS (Hess et al., 1997) algorithms, respectively. The equations of motion were solved using the Leapfrog scheme with a 1-fs time step. Subsequently, 20-ns constrained MD simulation was performed, using the dimensions and atom positions of the last snapshot of the NPT simulation. The settings for this run were kept the same as the previous section except that we turned the barostat off and the velocity scaling on.

For MD simulations of the RNA HJH, similar settings were used as for the riboswitch simulations. Here, we summarize the parameters that are different. For HJH, we used a PME set up with a grid spacing of 1.2 Å with an interpolation of order 4. The distance cutoff for non-bonded interactions and neighbor search was set to 11 Å. After equilibration, we performed 200 ns constrained MD at canonical ensemble during which heavy atoms of RNA were restrained, whereas water and ions were let free to move, and the output coordinate of this equilibration simulation was saved for the next step. To prepare the starting structures for the ensuing SAXS-driven simulations, we select conformations for each replica to reflect the structural diversity covering different regions of conformational space.

The theoretical SAXS profiles were computed from the simulation according to the theory of excess electron density (Park et al., 2009; Chen and Hub, 2014) (see Eqs.1-3 in Section 2). By following the methods described previously (He et al., 2021), we built a spatial envelope to enclose the ion and water near the biomolecule with a distance cutoff of d = 10 Å from the RNA surface (see Figure 1B) to ensure a bulk-like solvent at the envelope isosurface. As shown in Figure 1B, on the basis of the solute system (RNA + ion + water), the molecular envelope was constructed and extended to bulk. In addition, to accurately estimate the electron density of the bulk solvent and approximate the excluded volumes, we ran 20-ns-long MD simulations of buffer environment at the NVT ensemble. The buffer simulation contains the ion pairs and water molecules in a periodic box with dimensions and concentration matching the solute system. Similar to the solute systems, an envelope that extends to 10 Å was applied to the solvent system for buffer subtractions.

The agreement between experimentally and computationally acquired SAXS profiles was quantitatively evaluated using the linear χ ² metric χ 2 = 1 n − 1 ∑ i = 1 n I exp q i − I c o m q i σ q i 2 , (8) where n is the number of q points, and I _com (q _i ) and I _exp (q _i ) are the computed and experimental intensities at q _i , respectively. In Eq. 8, σ(q _i ) is the experimental error.

In this work, we demonstrated two possible ways to couple the experimental information with MD simulations—a direct coupling approach based on harmonic restraints and a maximum entropy–based approach. We used the first approach to analyze the SAXS profile of the SAM-I riboswitch that is assumed to be dominated by a unique conformational state. In comparison, the structural interpretation of the flexible HJH was accomplished using the maximum entropy approach.

During SAXS-driven MD (Chen and Hub, 2015), we used stochastic dynamics (SD) integrator to sample conformations. We removed the center of mass motion to keep the RNA atoms at the box center, while maintaining other settings from the previous section of GeneralMDsimulationsetup. The scattering amplitudes were averaged by giving more weight to recent conformations using an exponential decaying memory function (Chen and Hub, 2014). The memory time τ was set to picoseconds and enabled the convergence of buffer subtraction over the time interval covering the fluctuations of solvent and the RNA backbone. A time interval t of 20 ns was used for the switch function α(t) to avoid strong SAXS-derived bias and to allow sufficient convergence of I _com (q, t) before applying the SAXS-derived forces.

The computed SAXS curves were recorded every 2.5 ps for analysis. The uncertainty of solvent density was set to 0.5%. We used 1,500 q-vectors for orientational average. The structures were refined using the experimental data upto q = 2.94 nm ⁻¹, q = 3.15 nm ⁻¹, and q = 3.15 nm ⁻¹ for the models s1, s2, and s3, respectively. The number of restrained q points n _q was roughly determined according to the Shannon information theory (Shannon, 1949). We used n _q = 50 for our study on the SAM-I riboswitch structures. For the coupling potential, we chose k _c = 0.5. The convergence of simulations was evaluated frame by frame by χ ² metric (Eq. 8). The conformations giving good agreement with experiments (χ ² ≤ 2.0 in our case) are collected and used for further analysis.

Clustering of the data was done using the single linkage with a 0.1-nm cutoff, and the cluster centers were used to represent the corresponding ensembles. Root mean square fluctuation (rmsf) and the radius of gyration (Rg) of the RNA conformations, contact between the residues, and local Mg ²⁺ concentration were analyzed for conformations that show below χ ² threshold value. In the case of structure s3b, this value was set to 15, for all other structures to 10.

In parallel-replica SAXS-driven MD (Hermann and Hub, 2019) of the RNA HJH, all settings are the same as single-replica simulations except for certain molecule-dependent parameters. As aforementioned, a four-replica protocol was applied in the study of HJH. We used a time interval t of 10 ns for the switch function α(t), which allows sufficient convergence of I ̄ c o m ( q , R 1 , … , R N , t ) before introducing penalty. The structures were refined against the experimental data up to q = 2.9 nm ⁻¹, and the coupling strength k _c was set to 1.0 in this case. We used 800 q-vectors for the orientational average of real-time SAXS calculation. We analyzed the convergence of simulation frame by frame using χ ² metric and continued the simulation for at least another 100 ns upon reaching χ ² threshold ( ≤ 1.5). The conformations giving good agreement with experiments (χ ² ≤ 1.5) were selected and combined to build the conformational pool. The structures of the pool were clustered using the GROMOS (Daura et al., 1999) algorithm with a cutoff of 4.0 Å, and the cluster centers were used to represent each corresponding ensemble.

To describe the local Mg ²⁺ concentration, we selected the phosphate group of each residue. The concentration around each phosphate group is estimated by the following expression (Nguyen and Thirumalai, 2020): c ( X i ) = g X i r max c b u l k , (9) where X _i represents the phosphate group of residue i. g _Xi (r _max) is the maximum peak height of the normalized radial distribution function (RDF) between the cation and phosphate atom(s) on the RNA, and c _bulk is the bulk solution concentration.

To cross-validate the sampled conformations, we computed the average FRET efficiency and compared with the experimental data. The theoretical FRET efficiency is estimated on the basis of the work of Schuler et al. (2005): E F R E T = ∫ E r P r d r w i t h E r = 1 + r / R 0 6 − 1 , (10) where R₀ is set to be 60 Å after considering the linker length on the basis of previous works (Rothwell et al., 2003; Kalinin et al., 2010; Sindbert et al., 2011) and P(r) is the normalized distribution of distance (r) between the RNA surface where the FRET dyes attached in the experiments.

Using our methodology, we study the structural changes of metF-H2 SAM-I riboswitch from Thermoanaerobacter tengcongensis. In the presence of divalent ions and a ligand, this RNA folds into a four-way junction motif comprising a stem (P1; Figures 2A,B) and three stem-loop structures (P2–P4; Figures 2A,B). P1/P4 and P2/P3 form a co-axially stacked helix system where both helices are connected by a pseudoknot (PK) between P2 and P4, and joining regions J1/2 and J3/4. The linearity of P2 is interrupted by a kink-turn (KT), which facilitates the interactions in the PK region. As a result, further key tertiary contacts that connect P2 with J3/4 are formed. In addition, the crystal structure of the ligand bound complex shows a triple interaction between (A85–U64)•A24 and minor groove interactions between the P2 and A61-2 (Figure 2A). These network of interactions resulted in the three-dimensional structure shown in Figure 2B where the SAM ligand binds to a pocket that is situated between the helices P1 and P3 and the junction J1/2. The interaction of the adenine nucleobase of SAM with the asymmetric internal loop motif of the P3 helix is mainly characterized by hydrogen bonding, whereas the interaction with the minor groove of helix P1 is based on van der Waals surface complementarity. In addition, hydrogen bond formation can be found between the methionine moiety of SAM and J1/2.

Previous reported SAXS results by Stoddard et al. confirmed the similarity between the computed SAXS profile of the crystal structure of the ligand-bound SAM-I aptamer domain in the presence of Mg²⁺ with the solution SAXS profile (Stoddard et al., 2010). The profile of the solution structure was in good agreement with the theoretical SAXS curve derived from the crystal structure.

In addition to the ligand-bound structure, the SAXS study reports the aptamer structure in other important functional states, namely, the riboswitch without SAM in the presence of Mg²⁺ (apo-state, condition s2) and the apo-state without Mg²⁺ (condition s3). The reported Kratky plots of the scattering data show major differences depending on each of the conditions (Figure 2C). The changes in the SAXS profiles suggest structural differences in the riboswitch; however, an experiment-alone approach cannot resolve the atomic details.

To elucidate the implied structural differences reported in the SAXS profiles, we started with the SAXS profile of the SAM-bound riboswitch in the presence of Mg²⁺ (condition s1). Although the structure in this condition is already known, this exercise was aimed to establish our methodology and provide a benchmark against an independent data set. We have shown our results in Figure 3. Starting from the initial model adopted from the crystal structure, we computed the SAXS profile and monitor the error in SAXS as compared with the solution study. The χ ² value for the crystal structure was about 30 at the beginning, but this value rapidly dropped to χ ² ≈ 5 and started fluctuating around it, suggesting a convergence. We selected the lowest χ ² value (1.88) structure as the representative of the RNA at this condition (Figure 3B). To compare the structure derived from our methodology, with the initial model and the experimental result, we compared the SAXS profiles (Figures 3C,D). The intensity profile averaged from our pool generated after the sampling converged shows excellent agreement with the experimental profile, whereas the profile computed from the crystal structure shows deviations at high q vectors suggesting subtle structural changes in the RNA at short pairwise distances in the solution condition. To resolve the differences, we overlay the representative solution structure to the crystal structure (Figure 3E). Our methodology suggests a slightly more expanded structure for the RNA in solution with minor structural changes localized in the P1/P4 stack (Figure 3E). This result is consistent with the expectations, as crystallization conditions are known to lead the compact states for biomolecules. This result also verifies our methodology and lays the foundation for further investigations of the apo-states that we will discuss below.

Next, we show the structural ensemble of the apo-state of SAM-I riboswitch in the presence of Mg²⁺ (Figure 4). Time evolution of the error between experiment and simulation suggests a convergence after tens of nanoseconds. The lowest χ ² reported here is 1.12. Figure 4B displays the best fit model, whereas Figure 4C compares the average SAXS profile of the ensemble (blue) with the experimental data (gray), contrasted with the SAXS profile of the initial model (red).

The atomic model of SAM-I riboswitch in the apo-state derived from SAXS-driven MD is depicted in Figure 4B. Marked differences are evident when overlaid with the SAM-bound state (Figure 4D.) In the apo-state, we observe an overall expansion of the riboswitch dimensions. Interestingly, most of the changes are localized in the P1/P4 stack region. The changes in the P1/P4 stack likely result in the loss of the tertiary interaction between P1 and P3 domains that can alter (A85–U64)•A24 triple interaction (Figure 4A). Although the interaction between P1 and P3 is lost, both helices remain close to each other, suggesting that the initial collapse of the RNA is due to electrostatic forces and the binding of magnesium ions.

To analyze the riboswitch structure in the ligand bound and unbound states in more detail, we compute the contact map from each ensemble. In Figure 5, we show the occupancy of contact pairs in s1 and s2 conditions. The important domains of the riboswitch and residues that show marked differences in contacts are highlighted in the insets. Interestingly, despite the notable differences in the SAXS profiles (Figure 2C), the close contact pairs show similar patterns. The absence of the ligand results in losing or gaining contacts. For example, the contacts between the PK-region of P2 and J4/1 are lost in the absence of SAM (residues highlighted in cyan-blue and blue-yellow). Instead, the residue A85 forms a new contact with residue U26 in the s2 condition (residues highlighted in blue). In addition, new contacts between P1 and P3 are formed (residues highlighted in green). These new contacts likely narrow down the SAM-binding pocket surrounded by residues A6-C8 and A46-C47, and G56-G58 and G89. The structural details revealed by our model suggest that SAM binds to an already formed binding pocket.

Our data show that the riboswitch in the presence of Mg²⁺ ions is able to stabilize most of its tertiary contacts in the apo-state. To gain insights into how Mg²⁺ ions are stabilizing the s2 state, we computed local Mg²⁺ ion concentration around the oxygen atoms of the phosphate group of each residue. We compared cation distributions in the s1 and s2 states. Figure 6 shows the variation in the local Mg ²⁺ concentration as a function of residue identities. On the basis of this analysis, we identified two major binding sites for s1. These binding sites are located around the P2 region. The first one is around the residues G22-C30, whereas the second binding site is around the G35-C40 downstream of the KT. The major change between s1 and s2 occurs around the PK region. The loss of contacts around G22-C30 leads to the diminishing of PK-Mg²⁺ binding site in the s2 state. In addition to the major binding sites, we observe binding sites specific to either s1 or s2. For example, in s1, Mg²⁺ ions show strong localization around A70-A75 due to the increased charge density in P4-P1 contact. This binding site is not present in s2 as the contact is lost. In contrast, in the s2 state, we observe an s2 specific binding site around A53-G56 due to the newly formed tertiary contacts between P1 and P3 domains.

A recent experimental study by Sarkar et al. identified two similar Mg²⁺ binding sites, around U26 and A36 (Sarkar et al., 2021). In addition, their study shows a decrease in Mg²⁺ association at U26 below a certain SAM concentration confirming the observed changes of the KT–Mg²⁺ interaction site in our study. However, we could not verify the Mg²⁺ binding site reported between J1/2 and J3/4 in the s1 state (Stoddard et al., 2010; Sarkar et al., 2021) in our simulations. This is most likely due to the limitations in our approach. As the prime focus of this study is the RNA structure, we did not run our simulations long enough to observe convergence in ion distributions. The chelated Mg binding sites require long simulation times. Further study is underway to investigate this issue.

Next, we investigated the structural ensemble of the riboswitch when there is no ligand and no Mg²⁺ ions are present (condition s3). As physiological conditions always have Mg²⁺ ions, this exercise helps us to further understand the functional role of Mg²⁺ ions in modulating the active structure of the riboswitch. The experimental data in the absence of Mg²⁺ ions (Figure 2C) indicate an open structure for RNA at this condition. For that purpose, a partially folded RNA that retains its secondary fold is assumed to be a good starting point to model this state. To obtain a suitable starting structure, we gently disrupt the tertiary contacts either by pulling the two helical stacks apart or heating the structure at 350 K. The resulting conformations were clustered, and the theoretical SAXS profile for each unique structure was calculated to identify candidate structures that may give rise to an agreement with the experimental data. The details of getting initial models are described in the Methods (Section 2.3).

From this structural pool, we selected about 10 structures and sampled conformations from each initial model independently using SAXS-driven MD simulations. Four of these structures reached a χ ² value below a threshold (3a, χ ² = 5.46; 3b, χ ² = 12.35; 3c, χ ² = 7.08; 3d, χ ² = 7.39). The remaining simulations end up with high χ ² and hence were discarded. The structures proposed by SAXS-driven MD and their Kratky plots are shown in Figure 7 in comparison to the experimental profile and the initial model. This result highlights the advantage of our approach, namely, the capability to obtain a structural ensemble that represents the experimental SAXS data instead of providing only one possible structure. This is especially important for RNAs that possess structural flexibility and conformational heterogeneity.

The absence of the ligand and Mg ²⁺ ions is reflected in the overall size of the RNA. The average radius of gyration, R _g , measuring the overall size of a macroion shows a value of 2.14 ± 0.02 nm in the s1 state. This value increases to 2.39 ± 0.02 nm in the s2 state. Interestingly, the absence of binding partners leads to a major expansion of the structure, the R _g value becomes 2.75 ± 0.04 nm, in the s3 state. The resulting ensemble in the absence of binding partners is characterized by a larger gap between their P3 and P1 helices. The central part of the SAM-I riboswitch shows more open conformations. Similar to the s2 structure, here, the helix P4 is elongated. The similarities between s2 and s3 are likely due to a missing (A85–U64)•A24 triple interaction and the absence of short range interactions between helix P1 and P3 which happens in the absence of SAM. Despite the similarities between s2 and s3, there are additional changes unique to the s3 state; namely, an alternative base pairing of helix P4, the disruption of the base pairing between G82 and C69, and a 1-nt shorter loop structure due to an additional interaction between G74 and C78. As a result, P4 is less structured in the absence of Mg ²⁺ ions.

The structural changes due to binding partners are also reflected in the dynamics of the riboswitch. To elucidate the dynamical differences in each condition, we monitored the residue level RMS fluctuations (Figure 8) from our conformational pool. The conformations representing different clusters aligned to display the mode of fluctuations (Figure 8B). The results indicate that the fluctuations decrease with higher compaction of the RNA molecules, due to the presence of Mg ²⁺ ions and SAM. The region with the highest fluctuations was under all conditions the 5′ and 3′ end of the RNA situated in helix P1, followed by the P4 domain, and the loop of P3. The flexibility of RNA shifts from the loop to the stem region of the P4 helix in the absence of the ligand. Interestingly, the PK region remains rigid in all conditions, with some minor differences within the P2 helix.

Different from the riboswitch where the three-dimensional topology is relatively well restrained by tertiary contacts and ligands, the HJH motif (Figures 9A,B) possesses high structural heterogeneity. The flexible junction region allows its folding to be diverse, but at the same time, the base stacking of the junction region leads to distinct conformational sub-states. This high structural heterogeneity is modulated by solvent conditions, such as ion strengths and valence. Indeed, the conformations sampled in different salt conditions show detectable differences in the SAXS profiles for HJH RNAs (Chen et al., 2019; Chen and Pollack, 2019; Zettl et al., 2020; He et al., 2021).

Here, we investigate the conformational ensemble of HJH in varying ionic strengths of monovalent salt. In Figure 9C, we show the experimentally measured SAXS curves of RNA HJH represented by Kratky plots, I(q)q ² vs. q, at 50 mM (low), 100 mM (medium), and 500 mM (high) [KCl]. Typically, the Kratky plots highlight the scattering features at high q up to 0.3Å⁻¹. A more pronounced ridge in Kratky plots indicates the compactness of the molecule under study.

HJH’s sensitivity to salt concentrations is quite evident from the plots. Scattering profiles at low [KCl] and medium [KCl] show subtle deviations, especially at q ∈ (0.05 Å⁻¹, 0.1 Å⁻¹), which suggests conformational variations induced by increased [KCl]. Interestingly, such changes are amplified at high [KCl]. In particular, the reduction of the peak at mid-q regime [q ∈ (0.1 Å⁻¹, 0.15 Å⁻¹)] and the diminished peak at q ∈ (0.05 Å⁻¹, 0.15 Å⁻¹) indicates that the structure is even more extended at high-salt conditions. In other words, the raw SAXS data itself provide an interesting but counterintuitive observation—the cations are expected to soften electrostatic repulsion, which would lead to more compact states; however, the Kratky plot suggests more open conformations. The real-space projection of the SAXS profile could help to understand this observation.

By biasing simulations using the parallel-replica protocol (see details in Section 2.2), we obtained conformations that give rise to the experimentally consistent SAXS profiles. Here, we used four replicas as an optimum number to maximize the sampling and computational cost. Figure 10 summarizes the results for HJH simulations at different salt conditions. The structural search reached a plateau (χ ² ≤ 1.5) within 20 ns (Figures 10A–C) with computed profiles that align well with the experimental data (Figures 10D–F) with a χ ² value of 0.88, 0.85, and 1.35 for low, medium, and high [KCl], respectively. Interestingly, the χ ² scores are lower than the previous study (Chen et al., 2019) using the ensemble optimization method (EOM) (Bernadó et al., 2007). The success is likely due to the fact that, in SAXS-driven MD, the structures are committed to agreeing with the experimental signal; whereas in the re-weighting approaches, the structures are selected from a priori–generated pool. The lack of overlap between the experimental and computational pool may lead to artifacts in re-weighting approaches.

We then analyzed the trajectories of each salt condition that gives rise to good agreement with experiments. To visualize the structures, we align the bottom helix as shown in Figure 9B and project the HJH tip of the upper helix to a spherical coordinate system. The occupancy of conformations sampled is shown as a heat map (Figure 11). Interestingly, the conformations sampled by HJH show detectable differences depending on the ionic strength. In the low [KCl] (Figure 11A), we observed narrower distributions of helix orientations and less structural diversity reflected on the bending and twisting angles (see detailed definitions of orientation and bending in Figure 9B). At medium [KCl], the structural pool shifted to different bending angles retaining the twisting angle-similar to the low [KCl]. The conformations sampled at each salt condition are shown in (Figures 11D,E) for comparison. When it comes to high [KCl], however, the spherical heat maps reveal two distinct populations (Figure 11C). The HJH conformations now contain both extended (cyan) and bent states (blue) (Figure 11F). The bent states dominate the conformational pool at high-salt condition, yet the extended HJH states still contribute to the signal resulting in an average profile. Our visualization strategy resolves the conformational ensemble information hidden in the SAXS profiles. We conclude that the diminished peak at q ∈ (0.1 Å⁻¹, 0.15 Å⁻¹) of the SAXS curve is the result of this two co-existing populations.

To further assess whether the conformations sampled by this approach lead to over-fitting, one has to check the predictions using independent data sets. For that purpose, we computed the R _g and FRET efficiency for cross-validation. The measurements of R _g and E _FRET are available for the HJH system under study which helps to test the accuracy of the conformational ensembles generated by our approach. Figures 12A,B show our comparisons. Remarkably, the model shows the trends of R _g consistent with experiments. At low [KCl], both descriptors give rise to good agreement, which indicates the success of our strategy in capturing the conformational ensemble of HJH at low salt. However, at medium salt, the E _FRET values that we computed deviate from the measurement although the R _g derived from the computational model is in good agreement with the measurements. At high-salt conditions, our prediction of R _g is lower than the measurement, yet E _FRET value computed from the simulation is comparable to the experimental measurement. The disagreement between FRET and SAXS at medium salt could be attributed to the added complexity by chromophores. Our simulation does not include the chromophores explicitly, so their interactions with each other and RNA are ignored, which might explain the differences between simulation and experiment. In addition, our calculations are based on Eq. 10, which assumes the orientation of dye pairs to be isotropic and uncorrelated. This assumption may not be valid especially when the dyes come close to each other, which is the case in medium salt. Another possibility for the discrepancy might stem from the differences in the sensitivity of FRET and SAXS. For HJH geometry, FRET shows a significant change in efficiency when the two helices rotate, maintaining the same distance. However, SAXS that reports only pairwise distances may not distinguish the rotation. Our method that relies on SAXS may miss the rotation leading to the discrepancy. Further study is needed to resolve this issue.