Kinetics of Bulge Bases in Small RNAs and the Effect of Pressure on It

Due to their self-catalytic properties, small RNAs with bulge bases are hypothesized to be primordial molecules which could form elementary translation systems. Using molecular dynamics simulations, we study the binding propensity of small RNAs by calculating the free energy barrier corresponding to the looped out conformations of bulge bases, which presumably act as the binding sites for ligands in these small RNAs. We find that base flipping kinetics can proceed at atmospheric pressure but with a very small propensity. Furthermore, the free energy barrier associated with base flipping depends on the stacking with neighboring bases. Next, we studied the base flipping kinetics with pressure. We find that the free energy associated with base looping out increases monotonically as the pressure is increased. Furthermore, we calculate the mean first-passage time of conformational looping out of the bulge base using the diffusion of reaction coordinate associated with the base flipping on the underlying free energy surface. We find that the mean first-passage time associated with bulge looping out increases slowly upon increasing pressures up to atm but changes dramatically for atm. Finally, we discuss our results in the light of the role of hydration shell of water around RNA. Our results are relevant for the RNA world hypothesis.


Introduction
RNA molecules are very diverse both structurally and functionally [1,2]. Apart from having the regular helical purinepyrimidine base pairs, RNA molecules are also found to have many other secondary structures (motifs) such as loops, knots and bulges etc [3,4]. The presence of such structural motifs is found to play a role in binding of different molecules to RNA [5]. For many protein binding RNAs, it was found that the frequency of adenosine bulge at the binding site is very high. The presence of a bulge may change the conformational flexibility of an RNA [6,7] and hence more internal surface area of RNA is available for any chemistry. Moreover, the presence of the bulges does not only change the conformational flexibility but the bulges themselves may just flip out exposing the internal regions of a RNA to solvent and ligands. It has been shown that the bulge base looping out is highly sensitive to the bulge bases and their neighbors [6,7], which makes the question of generality of any picture of base bulge looping out difficult.
RNA is hypothesized as a primordial molecule due to its ability to form efficient catalysts and its similarity to DNA, where it can act as a molecular information machinery [8]. Yarus and coworkers have shown that RNAs, as small as 29 nucleotides, can undergo self-aminoacylation with aminoacyl adenylate (aa-AMP) as a substrate [9]. Lehmann and coworkers further studied the effect of base composition and length of the 39 extension on the aminoacylation rate of these small RNAs [10]. These studies indicate that these small RNAs may act like primitive tRNA.
Furthermore, due to rich chemistry of the reactive gases, and dissolved elements, a reminiscence of the early earth, hydrothermal vents are hypothesized to offer conditions viable for the origin of life [11,12]. Thermal vents naturally provide thermal gradients over small to very large length scales. Besides temperature gradient, the pressure near vents is much higher compared to atmospheric pressure. Hence, if a small primitive tRNA is selected at these thermodynamic conditions then they should also be able to carry out self-aminoacylation at high pressures.
The active site in the small tRNA-like molecules, that is believed to bind to aa-AMP, consists of AG/AG bulge in the stem [9,10]. The binding of aa-AMP proceeds by flipping of the bulge bases providing space for adenosine of the aa-AMP to stack in. Recent state of the art computer simulations of small RNAs have shed light on the bulge base looping out processes [13][14][15] in few of the possible cases. Specifically, these works have studied the free energy barriers associated with bulge base looping out process [13]. For example, A. Barthel and M. Zacharias studied the kinetics of bulge base looping out of a single uridine and adenosine bulge structures. Specifically, they calculated the free energy barrier corresponding to torsional deviation which measures the degree of looping out of bulge bases from the local helical plane. They find that the conformational free energy change in the case of adenosine bulge in a complete looping out process is higher by 1:5 kcal.mol 21 as compared to the uracil bulge [13], suggesting that in a base nonspecific binding process, a structure with single uracil bulge base would have higher propensity to flip out of helical plane. Although a wealth of literature is available on the base looping from the helical plane at ambient conditions, the changes in the kinetics of bulge base flipping is unexplored at conditions away from ambient conditions. Related is the question of the effect of RNA hydration on RNA kinetics. For binding to proceed, the bulge base has to flip out of the local helical plane, and so has to overcome both the bending rigidity and solvation energy barrier to solvate in water. The solvation of different substances in water is a widely studied problem [16][17][18][19][20][21]. Indeed, studies of apolar solutes in water shows an elliptic region in the pressure-temperature plane in which water behaves as a bad solvent and hence less solubility of these substances [22]. Moreover, It is known that water's hydrogen bond network and so the local structure of liquid water changes upon changing thermodynamic conditions, giving rise to anomalous changes in the dynamics and thermodynamics of water and aqueous systems. The structural stability and kinetics of proteins (where structural stability usually implies kinetically functional) is a function of pressure and temperature. The solvation barrier of apolar solutes in the case of proteins plays an important role both in the hydrophobic collapse of the polypeptides [17,22] as well as the stability of these structures as a function of pressure and temperature. The effect of the RNA hydration on the kinetics of RNA has not been given much attention.
In this paper, we first compare the propensity of bulge base looping out of a double strand RNA (dsRNA) for three different sequences: (i) with a single adenosine bulge in sequence 59GGGGAGG39-59CCCCCC39, and (ii) with a single adenosine bulge in sequence 59CCCCACC39-59GGGGGG39, and (iii) an AA-bulge in 59GGGGAGG39-59CCACCCC39. Next,we study the effect of pressure on the bulge base looping out for sequences (i) and (ii) with a single A-bulge. The organization of the paper is as follows. In section ''Methods'', we discuss the computational method, in section ''Results'', we present the results for the free energy barrier associated with torsional deviation of the adenosine bulge base from the local helical backbone at both atmospheric and elevated pressures. Then, we next present a mean first-passage time calculation associated with bulge base looping out process, and the effect of hydration on the kinetics. Finally we summarize and discuss our results in the ''Discussion'' section.

Methods
The energy minimized starting structures of three different RNAs (i) with a single A-bulge embedded between (GC) 2 (59-GGGGAGG-39/59-CCCCCC-39) ( Fig. 1(a)), (ii) with a single Abulge embedded between (CG) 2 ( Fig. 1(b)), and (iii) with an AAbulge (59GGGGAGG39-59CCCCACC39) embedded between (GC) 2 ( Fig. 1(c)) were created using the NAB/Nucgen module of the Amber10 program suite. The RNA structures were then solvated in 4000 TIP3P water molecules such that there was about 1:0 nm space left between the boundary of the box and the RNA molecule. For electro-neutrality, appropriate number of Na+ counter ions were added to the system. Energy minimizations were carried using the steepest descent (1000 steps) in GROMACS3.3.3 [23,24] with keeping the RNA atoms fixed. GROMACS is an all purpose molecular dynamics program which has been used to study various molecular systems [23]. The choice of GROMACS for our studies were mostly the convenince of it for doing dihedral restraint simultions. Any other programs that can provide a nucleic forcefield implementation, such as Amber [25] or NAMD [26] could have been used. Simulations were carried out using the GROMACS3.3.3 program with Amber99 force field. Equations of motion were integrated using a time step of 2 fs with periodic boundary conditions. The long range electrostatic interactions were treated with the particle-mesh-ewald (PME) method. After minimization, the system was slowly heated to temperature T = 300 K with positional restraints.
After the position restrained simulations, unrestrained Molecular Dynamics (MD) was carried out at four different values of the pressures P~1,1000,2000,3000 atm. Thermal equilibrium at a constant temperature T~300 K and different pressures was achieved using Berendsen thermostat and barostat respectively. The final equilibrated conformation was then used for as the starting conformation for umbrella sampling at different pressures.
To quantify the relative propensity of binding of RNA at different pressures, we chose the dihedral angle h (C19C19C19N1) (see Fig. 2) as the reaction coordinate for calculation of free energy. Since the conformational changes are very slow and an equilibrium sampling of torsional angles require much larger time scales than computationally feasible, we use biasing potential to calculate the free energy. Umbrella sampling method was used to calculate the relative free energy of bulge base looping out conformations of the RNA shown in Fig. 2. Harmonic umbrella biasing potentials U bias,hi~k (h{h i ) 2 with a force constant k~0:05 kcal.mol 21 .deg 22 were distributed uniformly along the reaction coordinate h at an interval Dh ref~5 o . Consecutive sampling windows of h were started from equilibrium structure of last run. For each values of the umbrella sampling window, we run a 2 ns simulation and record the value of h every 0:2 ps. The final potential of mean force (PMF) was calculated using the WHAM (weighted histogram method) [27]. The unbiased probability distribution P(h) at a given temperature T under WHAM is given by where N sim is the number of sampling window (simulations), n i is the number of counts in the bin associated with h, U bias,i is the biasing potential, and G i free energy from simulation i and is given by where h bins is the number of bins for the individual sampling window. The equations 1 and 2 are iterated to obtain the self consistent value of P(h). The value of P(h) depends on the time scale of simulations and hence long simulations are needed for a good convergence of the free energy.

Kinetics of bulge base at ambient and elevated pressures
In Fig. 3(a), we show DG(h) as a function of h for the RNA structure I with a single A-bulge at P~1 atm. Negative value of h corresponds to deviation towards minor groove while the positive values correspond to deviation towards major groove. We find that DG(h) has characteristic two minima centered around h&{20 o and h&30 o as reported in earlier studies of single A bulge [13]. Note that the definition of the torsional angle h is different from the one used Ref. [13] and hence different values of h. As we can see from Fig. 3 (a), the orientation of the A-bulge at more stable minimum is tilted slightly along the major groove while the second minimum at h&{20 o is presumably due to the base triplet formation with the neighboring bases. The free energy difference between these two minima is &3 kcal.mol 21 , suggesting that although h~30 o is relatively a more stable minimum configuration, the thermal fluctuations at T~300 K (&0:60 kcal.mol 21 ) is sufficient enough for the bulge to get displaced of the free energy minimum configuration. Due to a large free energy barrier (DG&5 kcal.mol 21 ) a complete looped out conformation of single A-bulge is less favorable and hence consistent with experimental fact that RNA with single A-bulge does not show appreciable aminoacylation (unpublished work). In Fig. 3 (b), we show DG(h) for a single A-bulge looping for structure II. In this case, the free energy barrier associated with complete looped state is smaller compared to structure I. Moreover, DG(h) exhibits two minima, one corresponding to the energy minimum stacked configuration and the other corresponding to the base-triplet formation at about 20 0 along the minor groove. In Fig. 3(c), we show DG(h) as a function of h for the RNA structure III with an AA-bulge at P~1atm, and T~300 K. Compared to single A-bulge in structure I (Fig. 3), we find that in the case of AA-bulge, DG for a complete flipped out configuration is smaller (&3 kcal.mol 21 ). Moreover, the secondary minimum as seen in the the case of single A-bulge (structure I, and II) is absent. Hence, the rate of base flipping would be enhanced for AA-bulge compared to single Abulge base in structure I.
We next studied the effect of pressure on the kinetics of bulge bases in structure I and structure II. Fig. 4(a) shows DG(h) for Abulge in structure I as a function of h for different P at T~300 K. We find that for pressures up to 2000 atm, F (h) has the characteristic two minima as we find in the case of P~1 atm. However, as the pressure in increased, the minimum at h&{20 o becomes shallower and disappears for Pw2000 atm, suggesting that at Pw2000 atm the base triplet formation of the bulge base with the neighboring bases does not occur during the looping process. Moreover, the free energy barrier between the two minima changes just a little upon increasing pressure for Pv2000 atm. For P~3000 atm, free energy barrier for the flipped out state changes drastically where DG&8:7 kcal.mol 21 , suggesting that the propensity of single A-bulge base looping out from the local helical plane would decrease upon increasing pressure and so the binding propensity of incoming ligands. Fig. 4(b) shows DG(h) for A-bulge in structure II as a function of h for different P at T~300 K. As the pressure is increased, the base triplet minimum seen at {20 0 disappears. The free energy barrier associated with a complete flipped out state of the bulge base progressively becomes larger with pressure exhibits sharp jump at P~2000 atm. At the largest pressure, P~300 atm, studied here, DG corresponding to complete flipped out state is &8 kcal.mol 21 . Comparing Fig. 4 (a) and (b), we find that the effect of pressure up to about 2000 atm is rather moderate on the kinetics although the rate of bulge base flipping decreases monotonically with pressure. Moreover, in both cases, whether the bulge base A is embedded with neighboring Guanine bases or Cytosine bases, the qualitative effect of the pressure on the kinetics remains the same.
In order to calculate the effective rate k eff of transition from a stacked to looped out conformation, we use Langevin equation [28,29]. Assuming the diffusion of the reaction coordinate h on an underlying free energy surface the dynamics of h is governed by where h is the torsional deviation and D is the diffusion constant, k B is the Boltzmann constant, and f (t) is the thermal noise with zero mean, Sf (t)T~0 and delta function correlation, Sf (t)f (0)T~2Dd(t). In the high friction limit, the probability r(h,t) of finding the system with reaction coordinate h after time t is given by the Smoluchowsky equation: where L is the Fokker-Planck operator given by L~L h e {bDG(h) DL h e bDG(h) and b~1=k B T. The mean first-passage time t(h i ) associated with crossing the barrier from any coordinate h i to final state h f is given by (see Fig. 5) where h r and h s denote the reflecting and absorbing boundaries respectively. We choose h f~{ 80 o as the final looped out state and initial state is chosen to the values of h where the free energy curve has the deepest minimum for respective pressures. The effective rate k eff of transition from I to III would then be given . The reflecting boundary was chosen to be at h~100 o where the relative free energy is &10k B T. Using Eq. 5, we calculate the value of t(h i ) for different pressures. We list the values of t(h i ) in table 1, where D is measured in deg 2 =sec. We find that t(h i ) increases upon increasing pressure within the error bars and increases sharply for P~3000 atm.

Hydration shell and base flipping of RNA
As we have seen in the sections above, the base flipping kinetics changes as the pressure is increased -namely, the free energy   Fig. 4(a) for structure II. The free energy barrier between the two minima slowly disappears upon increasing pressure, suggesting that the transient barrier that is produced at atmospheric pressures due to the partial base triplet hydrogen bonding of the bulge base with the neighboring bases is broken at higher pressures. Moreover, in both the cases DG corresponding to complete flipped out state shows a sharp change for Pw2000 atm. doi:10.1371/journal.pone.0042052.g004 barrier for the bulge base to flip out increases with pressure. Moreover, the transient barrier which presumably is due to the base triplet formation of the bulge base with the neighboring bases disappears at pressures Pw2000 atm. We note that Pw2000 atm is the pressure where most of the anomalies of liquid water disappears and also the pressure at which hydrophobic barriers for small molecules in water tend to vanish. Motivated by this we looked at the structure of the solvation shell (first hydration shell) of water around RNA for different pressures. We show a typical hydration shell around RNA in Fig. 6. The hydration shell is calculated by finding all the water molecules within a distance R C of RNA molecule. We choose R C~0 :223 nm as the first minimum in the radial distribution function of RNA and oxygen of water molecules (not shown here). We find that R C is independent of the pressure.
We first looked at whether the observed change in the pressure dependence of the kinetics is a result of ordering of water around RNA. To this effect, we calculated the OOO-bond angle w and its distribution P(w) of a central water molecule and its nearest neighbors in the hydration shell. For ordered liquid water w is very close to the tetrahedral angle 109:47 o . In Figs. 7 (a) and 7(b), we show P(w) for structures I, and II for pressures P~1,1000,2000,3000 atm, and T~300 K. For a comparison we also plot P(w) for bulk TIP3P water at P~1 atm and T~300 K. We find that, the second peak corresponding to more ordered water of the distribution P(w) shifts to smaller values of w, suggesting that the water monotonically disorders upon increasing pressure. Moreover, we do not find any significant sharp changes in P(w) which could be associated with sharp change in the free energy barrier observed at P~3000 atm.
To quantify the ordering/disordering of water molecules around RNA, we use the tetrahedral order parameterQ [30][31][32]. Tetrahedral order parameter Q quantifies how close a given water molecule and its first shell neighbors form a structure close to a tetrahedron. In general, Q k of k th molecules is defined as where, the indices i, and j run over all four neighboring molecules and w ikj the OOO-angle formed between the oxygens of central molecule k and neighbors i and j. Since, we only consider a thin hydration shell, the expression for ensemble average vQw can be written as vQw~1{ 9 4 In table 2, we list the average tetrahedral order parameter vQw of the hydration shell for different pressures. For a comparison, we also compute vQw for bulk water. Table 2 lists values of vQw for different pressures. We find that, vQw monotonically decreases upon increasing pressure and no sudden change in vQw is seen at P~3000 atm.
Since we did not see any sudden change in the ordering of hydration shell around RNA that might lead to the base flipping barrier observed at P~3000 atm, we next studied the hydration shell of the bulge. In Figs. 8(a) and (c), we show distribution P(n W ) of number of water molecules n W in the hydration shell of the bulge for structures I, and II respectively. Surprisingly, we find that at high pressures, the average number of water molecules in the first hydration shell of the bulge base increases, from an average of about 1:0 to 1:40 (see Fig. 8(b) and Fig. 8(d)). Moreover, we find that the distribution of water molecules around the bulge base shows significant probability of finding 3{4 water molecules. To this end, we suggest that the increased barrier of base flipping and the disappearance of the transient barrier in the free energy barrier is due to the presence of increased water molecules in the solvation shell around the bulge base. At high pressure, presence of more water molecules in the hydration shell of the bulge, which might have penetrated from the major groove, may have led to stabilizing the stacking of the bulge base. Increased free energy barrier at high pressures could be due to energetic factors such as hydrogen bonding with the increased water molecules in the hydration shell, or entropic factors. In general it is a contribution of both factors. In future, relative roles of enthalpic and entropic barriers must be investigated to shed more light on the sharp change in the kinetics observed at high pressure.

Discussion
In summary, in this paper we have investigated the effect of neighbor stacking and pressure on the kinetics of an Adenosine bulge base embedded between (GC) 2 in 59GGGGAGG39-59CCCCCC39 (structure I), and embedded between (CG) 2 in 59CCCCACC39-59GGGGGG39 (structure II), and AA-bulge embedded between (GC) 2 in 59GGGGAGG39-59CCACCCC39 (structure III). Specifically, we calculate the free energy barrier associated with the base looping out from the local helical plane of the RNAs. We find that a single A-bulge base embedded between (GC) 2 has a much larger free energy barrier of complete looping out compared to both AA-bulge embedded between (GC) 2 , and Abulge embedded between (CG) 2 . Among the three structures studied here, the free energy of looping out of A-bulge is minimum for A-bulge base embedded between (CG) 2 . Our results indicate the importance of stacking interactions with neighboring bases in determining the rate of base flipping. It further suggests that certain structural features of the bulge bases may facilitate the binding of small molecules to small RNAs. Usually, an A-bulge base with weaker neighbor stacking will have faster rate of flipping as compared to AA-bulge with the same neighbor stacking. We have studied only a few simple cases here. A complete underlying physical picture could only come from looking at various structural motifs and identifying the structural features that may increase or decrease the binding propensity of small RNAs.
Motivated by the fact that near the thermal vents the pressure is very high, we next studied the effect of pressure on the kinetics of A-bulge base for two configurations. We find that, upon increasing the pressure, the propensity or likelihood of base flipping decreases. At pressure Pw2000 atm, we see a sharp increase in the free energy barrier. Further, we calculate the time scale of flipping by mapping the problem of base flipping to a diffusion of Figure 7. Probability distribution function P(w) of the OOO-angle w in the first hydration shell at various pressures for (a) for structure I, and (b) for structure II. Note that a comparison with P(w) for bulk water at P~1 atm suggests that hydration shell is more disordered and this disorder increases monotonically upon increasing pressure. doi:10.1371/journal.pone.0042052.g007 reaction coordinate on an underlying free energy landscape from which we calculate the time scale of looping out of bulge base. We find that the time scale increases upon increasing pressure and changes dramatically at Pw2000 atm. We associate this behavior to increased solvation of the bulge base at high pressures. The effect of hydration, namely the increased hydration level of the bulge bases led to slowing down of the bulge flipping kinetics. It implies that the structural changes in water and the RNA hydration at high pressures may be relevant to self-catalytic properties of small RNAs. Indeed, water exhibits many anomalous behavior [33,34] as a function of pressure and temperature, including local structural changes in the liquid state [32]. It would be important to explore the kinetics of bulge bases as a function of pressure and temperature to see if certain range of pressure and temperature may facilitate the kinetics of bulge bases. Moreover, recent experimental studies suggest that in the presence of a temperature gradient, biopolymers such as RNA, DNA can be separated and selected depending on their shape, size and sequence apart from earlier known accumulation and depletion behaviors [35]. The vents create porous precipitates that have large connected pores. Adding to this is the fact that surfaces affect the local structure of water and hence may affect the hydration of RNAs leading to changes in the kinetics [36]. In future, it would be important to explore not only the effect of thermodynamic conditions but also the effect of surfaces on the kinetics.