System set up

All simulations reported hereafter were performed on two systems: the add aptamer domain complexed with adenine (Holo form) and without adenine bound (Apo form). In both cases we used the X-ray structure solved by Serganov et al. [PDB:1Y26] [7]. The ligand was removed to simulate the unbound state. MD simulations were performed using the Amber99 force field [30] refined with the parmbsc0 α/γ corrections [31] and the latest χ torsional parameters [24]. The general Amber force field [32] was used to parametrize the ligand. Partial atomic charges were assigned using the restricted electrostatic potential fit method [33] based on an electronic structure calculation at the HF/6-31G* level of theory performed with Gaussian03 [34]. The electrostatic interactions were calculated using the particle-mesh Ewald method [35] and bond-lengths were constrained with LINCS [36]. The systems were set-up following exactly the protocol described in Allnér et al. [20]: aptamers were solvated in a rhombic dodecahedron having 8 nm as box vector lenghts, with a Mg²⁺-H2O solution using approximately 11000 TIP3 molecules [37], and a recent parametrization for divalent cations [38]. The 5 crystallographic Mg²⁺ were initially kept at their respective position, whereas the additional 30 ions added to neutralize the system ([Mg²⁺] = 0.18 M) were randomly placed. A steepest descent minimization (150 steps) was performed followed by 200 ps of MD at constant temperature (298 K, using stochastic velocity rescaling [39]) and pressure (1 atm, using the Berendsen barostat [40]) with positional restraints on both RNA and ions so as to equilibrate water. This procedure was repeated first removing the constraints on the ions and then removing all the remaining constraints. Finally, 12 ns unrestrained simulations at constant volume were performed for each system.

Umbrella sampling

In order to compute the thermodynamic stability of the loop-loop interaction we employed US simulations with multiple harmonic restraints. The distance between the center of mass (CoM) of the backbone atoms of the two loops (L2: bases 20-26; L3: bases 48-54; Figure 2B) was used as a collective variable (CV). We will refer hereafter to this distance as L. 44 uniformly spaced reference values were taken in the range spanning from 12.5 to 34 Å, and restraints with stiffness k = 20 (kcal/mol)/Ų were employed. In the production phase of the US simulations each of the 44 windows was run for 5 ns. A very important issue in US simulation is the generation of the starting conformations. We here performed two independent US sampling simulations, using starting conformations generated with two different protocols (hereafter referred to as forward and backward). To generate the starting points for the forward US simulations we employed the same protocol used by Allnér et al. [20]. Namely, we ran a series of 44 short (0.25 ns) simulations with a stiffer restraint (k = 40 (kcal/mol)/Ų) keeping the CV at the 44 reference values, where each simulation was initialized from the last frame of the previous one. In this way, before each US window starting structure was sampled, we let the system equilibrate. The reference values were iterated allowing an increasing distance between the loops. The backward US simulation was initialized with an equivalent procedure but iterating the restraints in the opposite order, i.e. starting from the structure with undocked loops. In principle, if US simulations are converged, the result should be independent of the initialization procedure.

Analysis methods

The data were analyzed using the last 4 ns of each window. The potential of mean force (PMF) profiles were constructed using the weighted histogram analysis method (WHAM) [27] implemented by Grossfield [41] taking the CV values distribution resulting from the US simulations. This implementation of WHAM allows to compute errors with a bootstrapping procedure that assumes uncorrelated samples. To avoid artifacts due to possible correlations we instead adopted a blocking procedure. Namely, we split the final 4 ns of each trajectory in four blocks of 1 ns each and performed the WHAM calculation using only a single block from each simulations. The four resulting profiles are aligned at their CV starting value (12.1 Å for the forward profiles, 34.4 Å for the backward profiles) and error at each point is computed as the standard deviation among the four profiles divided by $\sqrt{4}$.

To define the number of stacking interactions and the number of base-pair contacts a local coordinate system was constructed in the centre of each six-membered rings, with the x axis pointing towards the C2 atom and the z axis orthogonal to the ring plane. The pairing and stacking relationship between two bases j and k is based on the vector r _jk , i.e. the position of ring center k relative to the coordinate system constructed on base j. The criteria for determining the canonical WC base pairs are: 1) the base pair must be AU or GC; 2) The relative position of the bases is compatible with the geometry of a WC interaction. The latter condition is considered satisfied when the product of the Gaussian function $\mathcal{N}\left(r_{jk};\mu ,\sigma \right)\times \mathcal{N}\left(r_{kj};\mu ,\sigma \right)>10^{-8}$. Mean µ and covariance σ were obtained from the empirical distribution of WC pairs in the crystal structure of the large ribosomal subunit [42]. The criteria for determining the non-canonical base pairs are: 1) the ellipsoidal distance ${\mathcal{D}}_{jk}\equiv \sqrt{x_{jk}^2/25+y_{jk}^2/25+z_{jk}^2/9}<\sqrt{2.5}$, and ${\mathcal{D}}_{kj}<\sqrt{2.5};2$) $\left|z_{jk}\right|$ and $\left|z_{kj}\right|<2\mathsf{\text{Å}};3$) it is not a WC pair. The criteria for determining the stacking base pairs are: 1) the ellipsoidal distance ${\mathcal{D}}_{jk}<\sqrt{2.5}$, and ${\mathcal{D}}_{kj}<\sqrt{2.5};2$) $\left|z_{jk}\right|$ and $\left|z_{kj}\right|<2\mathsf{\text{Å}};3$ ) $x_{jk}^2+y_{jk}^2<5\mathsf{\text{Å}}$ and $x_{kj}^2+y_{kj}^2<5\mathsf{\text{Å}}$. This procedure yields similar results compared to the MC-annotate software [43] and was shown to be useful for characterizing both structural and dynamical properties of RNA molecules [44]. The software used to perform this structural analysis is available online (http://github.com/srnas/barnaba).

Forward process

The analysis of the Holo forward and Apo forward US trajectories allowed the PMF for the disruption of the kissing complex to be computed. The resulting profiles are plotted in Figure 3 for both Holo and Apo systems. The PMF shows a minimum at $L\approx 12.5\mathsf{\text{Å}}$, corresponding to the initial structure. The free energy change upon disruption of the kissing complex for the Holo structure is $\mathrm{\Delta }G=52\pm 2$ kcal/mol. For the Apo structure the stability of the complex is largely reduced to $\mathrm{\Delta }G=35\pm 3$ kcal/mol. The stabilization of the kissing complex provided by the ligand can thus be estimated as $\mathrm{\Delta }\mathrm{\Delta }G=17\pm 3$ kcal/mol. To understand which are the interactions that are relevant for the kissing complex formation we analyze inter-loop pairings and inter-and intra-loop stacking interactions for each of the restrained simulations (Figure 4). For both the Holo and the Apo forms, at a $L\ge 16\AA$, only the two inter-loop WC base pairs (G25-C49, G26-C48) peculiar of the loop-loop interaction are still formed. On the contrary, all the non-canonical base pairs are disrupted. In the Apo structure the inter-loop WC pairings were irreversibly lost at $L>23\mathsf{\text{Å}}$, whereas in the Holo structure they are at least partially maintained until $L\approx 30\mathsf{\text{Å}}$. We also analyzed the rupture of stacking interactions, distinguishing intra-loop and inter-loop contacts. Inter-loop stacking behaves in a manner qualitatively similar to the inter-loop WC pairings, going to zero at a distance $L\approx 23\mathsf{\text{Å}}$ (Apo) and 30 Å (Holo). On the contrary, the intra-loop stacking interactions are still present when the kissing loop is disrupted, indicating that the internal structure of the two

loops is preserved during undocking. It can be observed that in the Holo simulation the number of intra-strand stacking slightly decreases (≈ 5) for $29\mathsf{\text{Å}}\underset{\backsim }{<}L\underset{\backsim }{<}30\mathsf{\text{Å}}$ because of the distortion in the structure induced by the one of the two inter-loop WC pairings. After this residual interaction is lost all the intra-loop stacking contacts are recovered.

Backward process

In order to better assess the convergence of the free energy landscape for kissing complex formation, we also reconstructed the PMF profiles of the Holo and Apo structure from US simulations initialized with the backward process (Figure 2). The forward and backward profiles were aligned at $L=34\mathsf{\text{Å}}$, since at that distance the starting structure of the backward process is equal to the final structure of the forward one. The free-energy change upon docking is estimated taking the difference between the minimum value of the PMF (Holo:$L\approx 19\AA$; Apo: $L\approx 16\AA$) and its value for the undocked structure ($L=34\AA$): for the Holo $\mathrm{\Delta }G=-3.2\pm 0.9$ kcal/mol, for the Apo $\mathrm{\Delta }G=-5.9\pm 0.6$ kcal/mol. Albeit negative, these numbers are too small and not compatible with the ones found in the forward process. This is a clear signature of hysteresis in the pulling procedure that strongly biases the initial starting points of the US simulation. The reason for this discrepancy can be better understood by performing a structural analysis of the interactions on the different US windows. As it can be seen in Figure 4, in the backward process the native WC base pairs are not reformed. In general, a few contacts are formed between the two loops but they are not enough to stabilize the kissing complex. To be sure that this is a systematic effect we also tried a few alternative settings for backward simulations. Results are presented in the Appendix.

Discussion and Conclusions

Our calculations provide quantitative and atomistic details on the mechanism of kissing loop breaking and formation in the add riboswitch aptamer domain. The results can be compared with those recently obtained by Allnér et al. [20] on the same system using the CHARMM36 force field [22, 23]. In particular the free-energy computed with the forward process has been obtained with an identical protocol so as to allow a fair comparison between the force fields. In our work the estimated stability of the kissing loop complex is $\mathrm{\Delta }G=52\pm 2$ kcal/mol (Holo) and $\mathrm{\Delta }G=35\pm 3$ kcal/mol (Apo), so that upon ligand binding $\mathrm{\Delta }\mathrm{\Delta }G=17$ kcal/mol > 0. On the contrary, Allnér et al. reported $\mathrm{\Delta }\mathrm{\Delta }G=-10$ kcal/mol < 0. The sign of ΔΔG indicates whether the ligand binding and the formation of the kissing loop complex are cooperative (positive) or anticooperative (negative). Results obtained with the two force fields thus interestingly lead to two opposite pictures.

Recent experiments probed the differential ligand affinity in aptamers with mutations hindering the formation of the kissing complex [14]. The change in affinity indicates a cooperativity between ligand binding and kissing complex formation. This stabilization has been estimated to be $\mathrm{\Delta }\mathrm{\Delta }G\approx 6$ kcal/mol. This number should be interpreted with caution since it is based on the assumption that the mutated aptamer mimics a ligand-bound state that is accessible to the wild type aptamer [14]. Results obtained with Amber force field are in qualitative agreement with this picture.

Recently, the thermodynamics of other stand-alone kissing complexes have also been studied using different biophysical techniques [45, 46]. In these two experimental works the stability of the loop-loop interactions were found to be in the range $8\underset{\backsim }{<}\mathrm{\Delta }G\underset{\backsim }{<}14$ kcal/mol. Stability depends on the exact sequence and set of intra-loop interactions, but is always on the order of ten kcal/mol. The estimated stability of the kissing loop complex in our calculation, namely $\mathrm{\Delta }G=52\pm 2$ kcal/mol and $\mathrm{\Delta }G=35\pm 3$ kcal/mol for Holo and Apo respectively, is thus much larger than expected. Results obtained with CHARMM indicate a lower ΔG for both the systems [20], in better agreement with experimental results, even if still overestimated. Our result could be affected by the known overestimation of stacking interactions in the Amber force field [47]. Additionally, we would like to point out that the overestimation found in our calculations could also be a consequence of difficult convergence in the US simulations. To test if the US simulation are effectively converged, we tried to recover the profiles from simulations initialized with the backward process, with a procedure inspired by two directional pulling in steered MD [48–50]. The discrepancy between forward and backward process is an index of high dependence of the PMF on the initialization procedure and poses some questions on the actual convergence of the US simulations. Similarly to steered MD, one can expect that simulations in the forward and backward process are respectively overestimating and underestimating the kissing complex stability. Optimal results can be obtained in steered MD by combining simulations performed with both protocols [49, 50]. We stress here that even if the forward simulations apparently recover the qualitative behavior of the general accepted model, they cannot be trusted for a quantitative estimation of the free-energy change. The fact that the backward process cannot reach the native docked state is a signature of a barrier in an orthogonal degree of freedom that is not properly sampled. A possible candidate is the barrier related to the desolvation of the loops, required to form the correct interstrand interactions. Additionally, we observe that pulling on the distance between the two loops does not necessarily induce the entropic reduction required upon docking. These issues are expected to affect both forward and backward pulling. Our simulation could not give an estimate of the additional barriers, but we can assume that these issues equally affects the Holo and the Apo systems. Thus, the converged ΔΔG upon ligand binding should be somewhere in between results from the forward and backward simulations. Thus we can expect the ΔΔG to be in a wide range between -2.7 kcal/mol and 17 kcal/mol, which is in qualitative agreement with already mentioned experiments [14].

The convergence problem is not related to the US method itself but to the difficulty of describing such a complex docking event using a single distance as a CV. This variable is not sufficient to drive the system through the appropriate transition states. This is likely due to the existence of additional barriers on hidden degrees of freedom (e.g. solvation). We believe that in order to reliably quantify the ΔΔG for this system with US or other biased sampling methods one should employ more complex CVs which are closer to the actual reaction coordinate.

In conclusion, in this work we addressed the formation of the kissing loop complex in the A-riboswitch aptamer by means of accurate molecular dynamics simulations in explicit solvent combined with enhanced sampling techniques in presence or absence of the cognate ligand. Results are compatible with experiments and suggest that the ligand stabilizes the kissing loop formation. However, our results also spot some weakness of the umbrella sampling method and call for calculations performed with more advanced techniques, which will be the subject of future investigations.