Mesoscopic modeling for nucleic acid chain dynamics

To gain a deeper insight into cellular processes such as transcription and translation, one needs to uncover the mechanisms controlling the configurational changes of nucleic acids. As a step toward this aim, we present here a mesoscopic-level computational model that provides a new window into nucleic acid dynamics. We model a single-stranded nucleic as a polymer chain whose monomers are the nucleosides. Each monomer comprises a bead representing the sugar molecule and a pin representing the base. The bead-pin complex can rotate about the backbone of the chain. We consider pairwise stacking and hydrogen-bonding interactions. We use a modified Monte Carlo dynamics that splits the dynamics into translational bead motion and rotational pin motion. By performing a number of tests, we first show that our model is physically sound. We then focus on a study of the kinetics of a DNA hairpin—a single-stranded molecule comprising two complementary segments joined by a noncomplementary loop—studied experimentally. We find that results from our simulations agree with experimental observations, demonstrating that our model is a suitable tool for the investigation of the hybridization of single strands.

Figure 1

(Color online) Mesoscale representation of the basic “units” comprising a nucleic acid chain: Phosphodiester bonds (green circles), sugar molecules (light blue pentagons), and nitrogenous bases (large colored circles). The diagram to the right illustrates the different units in our model: Sugar molecules (blue circles) are bonded by phosphates (green straight lines) to form the phosphate backbone of the nucleic acid (green box); colored pins represent the nitrogenous bases. Here and in the following figures, we use the following color coding: yellow stands for Thymine (T), purple stands for Guanine (G), orange stands for Adenine (A), and dark blue stands for Cytosine (C).

Figure 2

(Color online) Lattice imposed constraints for the nucleic acid chain configurations. (a) We use a three-dimensional triangular lattice. We constrain consecutive phosphodiester bonds to have an angle larger than or equal to 60° in order to mimic the stiffness of the sugar-phosphate backbone. The diagram illustrates the conformations allowed for two consecutive bonds in a three-dimensional triangular lattice. The black solid line and the black circles represent the reference bond and beads, respectively, and the purple dashed lines indicate the allowed conformations for the following bond. In the diagram, colored dots represent lattice sites which are nearest neighbors of the central blue dot. Different colors indicate the plane on which the site sits: top (red), middle (blue), and bottom (green). (b) The phosphodiester bond can be easily torqued, hence a sugar-base complex can take any spatial orientation provided it does not overlap with the phosphodiester bond. The diagram illustrates the ten possible orientations that a base (pin) can take (purple ellipses) for a given conformation of the polymer chain indicated by the black circles and black solid lines.

Figure 3

(Color online) Interaction energies (in arbitrary units) between pairs of nucleotides. We use the color code shown on the right to represent interaction strengths. We obtained the single-strand stacking energies from experimental data reported in Ref. 60 (Chap. 8), and the hydrogen bonding and cross-stacking energies from the duplex stacking enthalpies used in the Turner Rules [60] (Chap. 8).

Figure 4

(Color online) Chain motion on the lattice. In the panels, different colors indicate the different planes to which the sites belong: Blue for the central plane, red for the plane above, and green for the plane below. We label the central site “0” and we number the twelve neighboring sites from 1 to 12. Site “13” is an example of a next nearest neighbor of the central site “0” that is a nearest neighbor of sites “5” and “6.” (a) Initial configuration of a polymer chain comprising three monomers $(m_1,m_2,m_3)$ sitting on the nodes of a three-dimensional lattice. To illustrate our algorithm for the motion of the chain, we consider the motion of the monomer $m_2$. (b) Projection onto the central plane of the polymer configuration and the neighboring lattice sites. The color code is the same as in (a). Monomer $m_2$ can move with equal probability to any of the ten empty neighboring sites. (c) $m_2$ moves to site “7” (indicated by the black arrow). This site is a nearest neighbor of site “1” in which $m_3$ sits, but it is not a nearest neighbor of site “5” which is occupied by $m_1$. Since consecutive monomers in the polymer chain must occupy neighboring sites in the lattice, $m_1$ must “reptate,” i.e., it must move to a neighboring site which is also a neighbor of site “7.” These are sites {“0,” “6,” “11,” “13”} indicated by boxes. Purple boxes show sites which cannot be occupied {“0,” “6”}, because the final configuration would violate stiffness constraints. Black boxes indicate acceptable sites {“11,” “13”}. With equal probability, monomer $m_1$ can move to either of the two acceptable sites. The two possible final configurations are shown in panels (e) and (d).

Figure 5

(Color online) Single-strand DNA hairpins. (a) Sample configurations of a ssDNA hairpin, comprising 25 nucleotides with sequence $\mathrm{GCGTT}\text{−}{\mathrm{T}}_{15}\text{−}\mathrm{AACGC}$, on a three-dimensional triangular lattice. Spheres with different color indicate nucleotides with different bases: Orange for adenine (A), purple for guanine (G), yellow for thymine (T), and blue for cytosine (C). Note that the lattice symmetries are almost unnoticeable. (b) Schematic illustration of the transition between open and closed states for a hairpin loop. The hairpin switches between open/coil and closed/native states with characteristic rates $k_{\text{opening}}$ and $k_{\text{closing}}$.

Figure 6

(Color online) Sampling of the configuration space. Sampling of the configuration space at infinite temperature for polymers comprising (a) six and (b) eight monomers. Different color lines represent different values of the probability $p$. At high temperature, one expects all configurations to be sampled with equal rates (whose value is shown by the blue line). The average rate is the number of time steps in the simulation divided by the total number of configurations for the polymer. The number of configurations for a polymer with six (eight) monomers that sits on a three-dimensional lattice and satisfies the stiffness constraints indicated in Fig. 2 is 7500 (186 792). For each polymer size, we collected statistics for 5 000 0000 time steps. Our results demonstrate that the polymer samples conformation space more uniformly for smaller values of $p$. (c) Skewness $S\left(p\right)$ of the distribution of sampling rates of the conformation space of the polymer for different values of $p$. The skewness measures the asymmetry of the distribution. For perfect sampling, we expect the distribution to be normal, that is, $S=0$. For $p=0$, we find $S=0.35$ $(L=6)$ and $S=0.59$ $(L=8)$, in good agreement with this expectation. (d) Kurtosis $K\left(p\right)$ of the distribution of sampling rates of the conformation space of the polymer. The kurtosis measures the decay rate of the tails of the distribution. For a normal distribution, one has $K=3$. For $p=0$, we find $K=2.6$ for $L=6$ and $K=3.3$ for $L=8$, in good agreement with this expectation. Note that both $S$ and $K$ take smaller values for $p<0.01$ for the longer polymer. This suggests that as the length of the polymer increases, the differences in the distributions for small $p$ with respect to $p=0$ become smaller.

Figure 7

(Color online) (a) Average radius of gyration $⟨R_g⟩\sim L^{\nu }$ versus polymer length $L$. The solid line indicates the best fit to the expected power-law behavior $⟨R_g⟩\sim L^{\nu }$, obtaining $\nu =0.75\pm 0.07$. (b) Mean-squared displacement ${\mathrm{X}}_2\left(t\right)$ versus time for different polymer lengths $L=48–197$. Note that ${\mathrm{X}}_2$ grows linearly with time as expected in a diffusive process [78]. (c) By scaling the data in (b) by $D\sim L^{\alpha }$ with $\alpha =0.80\pm 0.13$, we are able to collapse all the data onto a single curve. The data displayed in the plots are averages over 5000 runs 100 000–150 000 time steps long using $p=0.05$.

Figure 8

(Color online) (a) Equilibrium energy of a hairpin with sequence A-TTTT-T. This hairpin has an internal loop comprising four T’s and a stem with a single base pair A-T. Differently from Fig. 3, we consider that the stacking energy of the T’s is zero. Under these conditions, a hairpin configuration can only take four energy values, $ϵ=-3$, when the hairpin is closed (A-T hydrogen bond is formed) and A is stacked with its neighboring T; $ϵ=-2$, when the A is stacked with its neighboring T and the hairpin is open; $ϵ=-1$, when the hairpin is closed but there is no stacking between the A and its neighboring T; $ϵ=0$, for all other cases. Under these conditions, the exact number of configurations for each energy level, $g_i$, can be computed. (b) Occupation number $⟨n_i\left(T\right)⟩$ of each energy level, $i=0$ to $-3$, as a function of temperature. Colored dots indicate the numerical results obtained from averages over 5 000 000 Monte Carlo steps (MCS) using the parallel tempering MC method [80]. Purple solid lines correspond to the theoretical expressions for the energy $e={\sum }_i{ϵ}_i{g}_i{e}^{-{ϵ}_i∕T}∕\mathcal{Z}$ in (a), and the occupation number $⟨n_i\left(T\right)⟩=g_i{e}^{-{ϵ}_i∕T}∕\mathcal{Z}$ in (b), where $\mathcal{Z}={\sum }_i{g}_i{e}^{-{ϵ}_i∕T}$ is the partition function. Note the excellent agreement between theoretical predictions and simulation results.

Figure 9

(Color online) Test of the equilibrium properties of a hairpin with sequence GGATAA-${\mathrm{T}}_4$-TTATCC. We performed simulations using the parallel tempering method [80] for the cases $p=0.01$, 0.05, and 0.1. Results correspond to averages over 10 000 000 MCS. Simulation temperatures in the range [0,1] were mapped into absolute temperatures using the conversion factor $T_m^{\mathit{Sim}}∕T_m^{\mathit{MFOLD}}=1.71\times {10}^{-3}$ obtained in Fig. 10 for the ionic conditions $\left[{\mathrm{Na}}^{+}\right]=1\mathrm{M}$ and $\left[{\mathrm{Mg}}^{++}\right]=0\mathrm{M}$. (a) Specific heat as a function of temperature calculated as (i) the derivative of the energy with respect to temperature $dE∕dT$ (solid lines), and (ii) $c=(⟨E^2⟩-{⟨E⟩}^2)∕T^2$ (symbols) [77]. At equilibrium, fluctuation-dissipation relations must be fulfilled and the two methods must lead to equal estimates. This is indeed what we observe. Furthermore, note that the agreement between the two methods is excellent even in the melting region when the heat capacity has its peak. We also checked that at low temperature the hairpin reaches its minimum energy configuration. Thus, this test demonstrates that we reach equilibrium in our simulations and that the equilibrium properties that we measure are correct. (b) Melting curve for the values of $p$ considered in (a). Note that the curves are insensitive to the specific value of $p$ in the range considered. As expected, the fraction of broken bonds in the stem goes from zero at low temperatures (where the low energy of the closed/native state dominates the partition function) to 1 at high temperatures (where entropy dominates). Blue dashed lines indicate the melting temperature in both plots: At the specific heat peak and at the point where $\theta =0.5$. We obtain for the two cases $T_m=341\pm 3\mathrm{K}$.

Figure 10

(Color online) Interaction Testing: Comparison of results from simulations with theoretical values. (a) Ratio of the melting temperature estimates, $r\equiv T_m^{\mathit{Sim}}∕T_m^{\mathit{MFOLD}}$, obtained with our model and with Zuker’s DNA folding server [1, 81] for different hairpins with stems of four, five, and six bases long. For all hairpins, the loop comprises four T’s. Melting temperatures from the folding server correspond to the following ionic concentrations: $\left[{\mathrm{Na}}^{+}\right]=1\mathrm{M}$ and $\left[{\mathrm{Mg}}^{++}\right]=0\mathrm{M}$. Note that in the interaction tables of Fig. 3, we do not consider any difference between oligonucleotides starting at the $5^{\prime }$ or the $3^{\prime }$ end. The reason for this modeling choice is that the fluctuations observed in the melting temperature obtained in each case, $5^{\prime }\text{to}{3}^{\prime }\left[T_m\left(5^{\prime }\right)\right]$ and $3^{\prime }\text{to}{5}^{\prime }\left[T_m\left(3^{\prime }\right)\right]$ sequences, were considerably smaller than the fluctuations of the whole data set. Therefore, the melting temperature used in the data shown corresponds to the average of both $T_m$’s. The red solid lines correspond to the average factors for hairpins with stems comprising four, five, and six base pairs. Note that the fluctuations are quite small. (b) Relative fluctuations $ϵ=(r-\overline{r})∕r$ of the ratio of temperatures with respect to the mean $(\overline{r}=1.71\times {10}^{-3})$. The red line represents the mean $\overline{ϵ}=-0.025$ and the gray band represents the region within one standard deviation of the mean. (c) Normalized distribution of $ϵ$. The red line corresponds to a Gaussian fit of zero mean and standard deviation $\sigma =0.15$. We obtained the melting temperatures from parallel tempering Monte Carlo simulations for temperatures in the range 0.06 to 1, and performing averages over $3\times {10}^6$, $5\times {10}^6$, and $10\times {10}^6$ MCS for hairpins with stems comprising four, five, and six base pairs, respectively. All data correspond to the case $p=0.05$, but we found no significant changes for different values of $p$. The analysis of different ionic conditions yields similar fluctuations but different conversion factors, suggesting that salt concentration and temperature play a similar role in our model.

Figure 11

(Color online) Kinetic properties. (a) Relaxation rates $k_r$ for the hairpin GGATAA-${\mathrm{T}}_4$-TTATCC versus inverse temperature. Black circles correspond to the data obtained from the simulations while red squares correspond to experimental data obtained by Ansari et al. with the same hairpin for two different experimental setups [16, 18]. The dashed line corresponds to the fit to a two-state model with Arrhenius dependence of the relaxation rates; cf. (b). Our data were obtained by averaging over different dynamical histories $(\sim 200)$ for temperatures in the range $T_m\pm 0.07T_m$. We run simulations with a wide range of MCS (1 500 000 to 1 800 000) to make sure that our estimates of the relaxation rates converged. Simulation temperatures have been rescaled by a factor $1.82\times {10}^{-3}$ to convert the temperature into Kelvin. Note that we use a different factor from that obtained in Fig. 10 to adjust to the experimantal conditions used in Ref. 16. In order to convert simulation rates to experimental rates, one needs to use a factor of approximately ${10}^{-9}$. This value suggests that the time scale of a single MC step is approximately one nanosecond. (b) Two-state analysis. $k_r∕[1+K_{\mathit{eq}}^{-1}\left(T\right)]$ (black dots) versus inverse temperature. $K_{\mathit{eq}}$ is known from equilibrium measurements (Fig. 9) as $K_{\mathit{eq}}\left(T\right)=1∕\theta \left(T\right)-1$. The solid red line corrsponds to the fit to the expression in Eq. (4). The two fitting parameters are $E_a$, the activation energy, and $A$, a phenomenological constant rate. The parameters for the best fit are $E_a=-3.5\pm 0.4\mathrm{kcal}∕\mathrm{mol}$ and $A=30\pm 2{\mathrm{s}}^{-1}$.

Figure 12

(Color online) Comparison of two different Monte Carlo schemes for the hairpin of sequence GGATAA-${\mathrm{T}}_4$-TTATCC. (a) Melting curve. Fraction of broken bonds $\theta$ with respect to temperature obtained from parallel tempering simulations averaging over $2\times {10}^7$ MCS. Black symbols correspond to the Monte Carlo scheme defined for the model (Table ) and red symbols correspond to results for scheme B (Table ). Note that the curves are practically undistinguishable. (b) Energy versus temperature. Black and red symbols correspond to parallel tempering results, following the color code in (a). Green and blue symbols correspond to single temperature simulations in the transition region. Each point corresponds to averages over 200 different histories 1 200 000 MCS long. Green symbols correspond to the scheme defined for our model and blue symbols correspond to scheme B. Note that while results for the model scheme show an excellent agreement with the parallel tempering results, results for scheme B show average energies far below the parallel tempering curve.

Figure 13

Interaction optimization. Comparison of results from simulations with theoretical values for different sets of interactions. Ratio of the melting temperature estimates, $T_m^{\mathit{Sim}}∕T_m^{\mathit{MFOLD}}$, obtained with our model and with Zuker’s ${\mathrm{T}}_m$ server [1, 81] for different hairpins with a stem comprising four base pairs and a loop comprising four T’s. Different symbols correspond to changing some of the interaction strength values shown in Fig. 3 to which black symbols (엯) correspond. Blue symbols (▵) correspond to reducing the $\mathrm{G}-\mathrm{C}$ bond strength. Yellow symbols (◁) correspond to increasing the $\mathrm{A}-\mathrm{T}$ bond strength. Green symbols (◇) correspond to decreasing the GC stacking strength. Red symbols (◻) correspond to performing the three previous changes simultaneously. Note these changes of the parameters result in a decrease of the fluctuations, the minimum variance corresponding to changing the strength of AT bonds, GC bonds, and GC stacking interactions at once. We obtained the $T_m$’s from parallel tempering [80] simulations averaging over $5\times {10}^6$ MCS, with $p=0.05$.