The Thermodynamic and Kinetic Properties of the dA-rU DNA-RNA Hybrid Base Pair Investigated via Molecular Dynamics Simulations

Abstract

DNA-RNA hybrid duplexes play essential roles during the reverse transcription of RNA viruses and DNA replication. The opening and conformation changes of individual base pairs are critical to their biological functions. However, the microscopic mechanisms governing base pair closing and opening at the atomic level remain poorly understood. In this study, we investigated the thermodynamic and kinetic parameters of the dA-rU base pair in a DNA-RNA hybrid duplex using 4 μs all-atom molecular dynamics (MD) simulations at different temperatures. Our results showed that the thermodynamic parameters of the dA-rU base pair aligned with the predictions of the nearest-neighbor model and were close to those of the AU base pair in RNA. The temperature dependence of the average lifetimes of both the open and the closed states, as well as the transition path times, were obtained. The free-energy barrier for a single base pair opening and closing arises from an increase in enthalpy due to the disruption of the base-stacking interactions and hydrogen bonding, along with an entropy loss attributed to the accompanying restrictions, such as torsional angle constraints and solvent viscosity.

1. Introduction

DNA-RNA hybrids, consisting of a DNA strand and an RNA strand, are important intermediates in numerous biological processes such as Okazaki fragment formation during DNA replication and the reverse transcription process [1,2]. Around 59% of human genes contain sequences capable of forming such hybrids [3]. These hybrids play major roles in DNA replication, transcription, methylation, repair, telomerase replication, and reverse transcription [1,2,4,5,6,7,8,9,10], all of which rely on the formation and cleavage of these molecules. Recognized by RNase H, these hybrids can be specifically hydrolyzed on the RNA strand without affecting the complementary DNA strand [11,12]. Furthermore, several human diseases, including certain cancers and neurological disorders, are linked to the formation of DNA-RNA hybrids [13,14,15,16,17]. More recently, the application of gene editing techniques like CRISPR-Cas9 has sparked renewed interest in understanding DNA-RNA hybridization [18,19,20]. The biological functions of DNA-RNA hybrid duplexes are closely associated with their mechanical properties, particularly their kinetics and conformational changes. In the process of structural transitions, the formation and disruption of individual base pairs are key steps in many biological processes. Therefore, investigating the thermodynamic and kinetic properties of individual base pair opening and closing is crucial for a quantitative understanding of the biological functions of these hybrids.

The opening and flipping of the base pairs in RNA and DNA duplexes have been studied extensively using various experimental techniques, including NMR spectroscopy [21,22,23,24,25,26], fluorescence correlation spectroscopy (FCS) [27,28], and single-molecule forces [29]. However, limitations in these techniques leave some fundamental issues, such as the dynamics of single base pair opening and base flipping, unresolved. Therefore, all-atom molecular dynamics (MD) simulations have been widely employed to investigate the base pair opening in RNA and DNA duplexes [30,31,32,33,34,35,36,37,38,39,40,41,42,43,44]. For example, Pan et al. used MD simulations to study structural fluctuations in RNA and DNA duplexes, discovering that RNA base pairs opened into the major groove, whereas DNA showed no significant tendency for base pair opening [36]. Similarly, Zgarbova et al. explored base pair fraying in both DNA and RNA, observing that AT/AU terminal base pairs were significantly less stable than GC terminal pairs [43]. Quantitative studies have further explored the average lifetimes, barriers, and transition rates for forming and opening base pairs in both RNA and DNA duplexes [35,40]. In contrast, there has been less quantitative research, both experimentally and theoretically, on base pair opening in DNA-RNA hybrid duplexes. Experimental studies have demonstrated that the stability of hybrid base pairs differs from that of RNA and DNA base pairs. For instance, with the same base sequence in nucleic acid double helices, the stability of the rA-rU base pair is similar to that of the corresponding rA-dT base pair but differs from the stabilities of the dA-dT and dA-rU base pairs in the other two duplexes. Specifically, the dA-rU base pair is less stable than the rA-dT base pair, following the order: dA-rU < rA-rU < dA-dT [45,46,47]. The thermodynamic parameters of hybrid base pairs were first experimentally obtained by measuring the melting temperatures of hybrid duplexes [48], but these parameters were later revised due to their dependence on sequence and ion concentration [49]. Additionally, some fundamental problems, including the barrier of opening and forming the base pair and the quantitative results on the transition rates during base pair formation, remain unresolved.

In this study, using all-atom MD simulation, we quantitatively characterized the thermodynamic and kinetic parameters of the dA-rU base pair in DNA-RNA hybrid duplexes by directly simulating the opening–closing switch of the base pair. The computed results showed that the thermodynamic parameters of DNA-RNA hybrid duplex base pairs were consistent with those of the nearest-neighbor model [49] and close to the thermodynamic parameters of the AU base pair in RNA [40]. Additionally, the temperature dependence of the transition rates for the base pair’s opening and closing processes was obtained.

2. Results and Discussion

2.1. Definition of Closed State (cs), Open State (os), and Transition State (ts)

As in our previous analyses [35,40,41,42], the terminal base pair dA-rU would undergo the closing–opening transition process (see Figure 1) through an intermediate transition state. The corresponding conformations can be categorized into three states: open, closed, and transition. These states are identified based on the time-dependent root mean square deviation (RMSD) for the terminal bases dA and rU, relative to their initial structure, as well as the backbone torsion angle ζ. This torsion angle refers to the dihedral angle formed by the four atoms C3′(i)-O3′(i)-P(i+1)-O5′(i+1), where i represents the ith nucleotide in a polynucleotide chain, as shown in Figure 2.

For the closed state, the RMSD value is centered around 0.7 Å and the torsion angle ζ is approximately −75° (−50° to −100°), as shown in Figure 2. In this conformation, the two terminal bases, dA and rU, exhibit slight vibrations around their initial positions, with the base pairing and the base-stacking interactions with neighboring nucleotides remaining largely intact. In the open state, the RMSD value varies from 2 to 13.5 Å, and the torsion angle ζ centers around 50°, as shown in Figure 2. In this conformation, the two terminal bases, dA and rU, shift away from their initial positions and flip into the solvent, leading to the disruption of base-pairing and base-stacking interactions. For the transition state, the RMSD value exceeds 2 Å and has a very short residence time, although the torsion angle ζ remains within the closed state region. Simultaneously, the terminal bases flip outward into the solvent. Depending on the transition pathway, the transition state can be further classified into two types, “ctc” (transitioning from the closed state and returning to the closed state) and “oto” (transitioning from the open state and returning to the open state), as shown in Figure 3.

2.2. Thermodynamic Properties of the Terminal dA-rU Base Pair

At each simulation temperature, the population distributions of the closed, open, and transition states were calculated using the Formula (1)

$$p_{cl}={\tau }_{cl}/\tau =\sum _{i=1}^{N_{cl}}{\tau }_i^{cl}/\tau , { p}_{op}={\tau }_{op}/\tau =\sum _{i=1}^{N_{op}}{\tau }_i^{op}/\tau , p_{ts}={\tau }_{ts}/\tau =\sum _{i=1}^{N_{ts}}{\tau }_i^{ts}/\tau$$

(1)

where $\tau$ is the total simulation time and ${\tau }_{cl}$, ${\tau }_{op},$ and ${\tau }_{ts}$ represent the total time spent in the closed, open, and transition states, respectively. $N_{cl}$, $N_{op}$, and $N_{ts}$ denote the total number of snapshots in which the corresponding conformations reside in the closed, open, and transition states. ${\tau }_i^{cl}$, ${\tau }_i^{op}$, and ${\tau }_i^{ts}$ refer to the ith duration the conformations remain in the respective state. The probability of occupying the closed state at different temperatures during the simulation is shown in Figure 4a. Each point represents the probability calculated over time intervals from the beginning up to the respective simulation time. When the simulation time exceeds 3500 ns, the occupied probability of the closed state has a stable value. To ensure that the system reaches equilibration between the open and closed states, the simulation time is set to 4000 ns for each temperature. The probabilities of the closed, open, and transition states at the end of simulations are provided in

Table 1. The average lifetime τ_ave (ns) and the total number N of occurrences of conformations at closed, open, and transition states at different temperatures (K) during the 4000 ns simulation.

Temperature	Closed State (cs)			Open State (os)			Transition State (ctc)			Transition State (oto)
T(K)	${\mathit{\tau }}_{\mathbf{a}\mathbf{v}\mathbf{e}}^{\mathbf{c}\mathbf{l}}$(ns)	Ncl	pcl	${\mathit{\tau }}_{\mathbf{a}\mathbf{v}\mathbf{e}}^{\mathbf{o}\mathbf{p}}$(ns)	Nop	Pop	${\mathit{\tau }}_{\mathbf{a}\mathbf{v}\mathbf{e}}^{\mathbf{c}\mathbf{t}\mathbf{c}}$(ns)	Nctc	Pctc	${\mathit{\tau }}_{\mathbf{a}\mathbf{v}\mathbf{e}}^{\mathbf{o}\mathbf{t}\mathbf{o}}$(ns)	Noto	Poto
370	27.41	78	0.53	24.18	77	0.47	0.23	435	0.025	0.205	663	0.034
380	21.21	89	0.47	23.74	89	0.53	0.225	441	0.025	0.201	736	0.037
390	16.16	102	0.41	23.06	102	0.59	0.216	447	0.024	0.195	762	0.037
400	12.95	114	0.37	22.14	114	0.63	0.211	462	0.023	0.188	821	0.039

.

Since the probabilities of the closed and open states are much higher than those of the transition states, the system was initially modeled using a simple opening–closing two-state model. According to the equilibrium probability distribution, the probability of each conformation among all base pairs in the hybrid can be calculated using the equation $p_i=\exp \left(-G_i/k_{B}T\right)/\sum _j\exp \left(-G_j/k_{B}T\right)$, where $G_i$ is the free energy of the ith conformation when the system reaches equilibrium, $k_B$ is the Boltzmann constant, T is the absolute temperature, and $\sum _j\exp \left(-G_j/k_{B}T\right)$ is the partition function. The free-energy differences $(\Delta G)$ between the two states can then be determined as (2)

$$\Delta \mathrm{G}=-k_{B}Tln({\displaystyle \frac{p_{op}}{p_{cl}}})$$

(2)

where p_op and p_cl represent the probabilities of the open and closed states, respectively; $k_B$ is the Boltzmann constant, and T is the absolute temperature. As shown in Figure 4b, the free-energy differences $(\Delta G$) exhibit a linear relationship with the reciprocal of temperature $\left(1/T\right)$. Using the Gibbs free-energy equation $\Delta G=\Delta H-T\Delta S$, where $\Delta H$ and $\Delta S$ are the enthalpy and entropy changes during the transition between the two states, the thermodynamic parameters of the terminal dA-rU base pair in hybrid duplexes were determined: $\Delta H=-6.83 \mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{m}\mathrm{o}\mathrm{l}$ and $∆S=-18.2 \mathrm{e}\mathrm{u}$. These values are consistent with those predicted by the nearest-neighbor model [49] and are close to the thermodynamic parameters ($\Delta H=-7.3 \mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{m}\mathrm{o}\mathrm{l}$, $∆S=-18.5 \mathrm{e}\mathrm{u}$) [40] for the AU base pair in RNA. These findings also align well with experimental data on thermal stability [45,46].

2.3. The Kinetic Mechanism of the Terminal dA-rU Base Pair

The average lifetimes of the closed, open, and transition states can be calculated through ${\tau }_{ave}={\sum }_{i=1}^N{\tau }_i/N$, where ${\tau }_{ave}$ is the average lifetime, ${\tau }_i$ is the ith lifetime of the conformation in its ith occurrence, and N is the total number of occurrences of the base pair in each state. The average lifetimes of the three states are listed in Table 1. It is clear that the average lifetime of the transition state is significantly shorter than that of the closed and open states. In the context of the opening–closing two-state model, the opening rate $k_{-}$ (from the closed state to the open state) and the closing rate $k_{+}$ (from the open state to the closed state) can be obtained using the formulas $k_{-}=1/{\tau }_{ave}^{cl}$, $k_{+}=1/{\tau }_{ave}^{op}$, where ${\tau }_{ave}^{cl}$ is the average lifetime of the closed state and ${\tau }_{ave}^{op}$ is the average lifetime of the open state. Figure 5a illustrates the temperature dependence of the average lifetime of the closed and open states. It can be observed that the average lifetime of the closed state exhibits a strong dependence on temperature, while that of the open state shows only a weak temperature sensitivity. This finding is consistent with both theoretical results [35,40,41,42] and the experimental findings [50] that the folding and unfolding rates have different sensitivities to temperature.

Under the three-state model (closed, open, and transition states), the transition rate from the transition state to the open state, and $k_{t\to o}$, is equal to the reciprocal of the average lifetime of the transition state oto: $k_{t\to o}=1/{\tau }_{ave}^{oto}$. Similarly, the transition rate from the transition state to the closed state, $k_{t\to c}$, is given by $k_{t\to c}=1/{\tau }_{ave}^{ctc}$. The transition path time from the open state to the closed state, $t_{tp}^{o\to c}$, and from the closed state to the open state, $t_{tp}^{c\to o}$, can be expressed as $t_{tp}^{o\to c}={\tau }_{ave}^{ctc}$ and $t_{tp}^{c\to o}={\tau }_{ave}^{oto}$. During the transition between the closed and open states of the base pair, the transition path times are much shorter than the lifetimes of the closed and open states and show only weak temperature-dependent behavior (see Table 1). Based on the transition-state theory [51,52,53,54,55,56], the average residence time in the open (t_op) (3) state, closed (t_cl) (4) state, and the transition path time $t_{tp}^{o\to c}$ (5) and $t_{tp}^{c\to o} (6)$, can be calculated as follows:

$$t_{op}={{\displaystyle \frac{1}{k}}}_{+}={\displaystyle \frac{2\pi }{\beta D^{*}{\omega }^{*}{\omega }_o}}\exp \left(\beta \Delta G_o\right)$$

(3)

$$t_{cl}={{\displaystyle \frac{1}{k}}}_{-}={\displaystyle \frac{2\pi }{\beta D^{*}{\omega }^{*}{\omega }_{\mathrm{c}}}}\exp \left(\beta \Delta G_{\mathrm{c}}\right)$$

(4)

$$t_{tp}^{o\to c}={\displaystyle \frac{1}{\beta D^{*}({\mathsf{\omega }}^{*}{)}^2}}\ln \left(2e^{\gamma }\beta ∆G_o\right)$$

(5)

$$t_{tp}^{c\to o}={\displaystyle \frac{1}{\beta D^{*}({\mathsf{\omega }}^{*}{)}^2}}\ln \left(2e^{\gamma }\beta ∆G_c\right)$$

(6)

where $\beta =1/k_{B}T$, with $k_B$ being the Boltzmann constant and T the absolute temperature. $D^{*}$ represents the diffusion coefficient at the top of the free-energy barrier; ${\left({\omega }_o\right)}^2$, ${\left({\omega }_c\right)}^2$, and ${\left({\omega }^{*}\right)}^2$ refer to the curvatures of the free-energy surface in the open state, the closed state, and at the barrier, respectively. $\gamma$ is Euler’s constant, $∆G_o$ represents the free-energy barrier height from the open state to the closed state, and $∆G_c$ represents the free-energy barrier height from the closed state to the open state. Equations (3) and (4) are attributed to Kramers [53], while Equations (5) and (6) are derived from Szabo [54,55]. Szabo’s equation makes the same assumptions and approximations as Kramers’ with regard to the underlying physics. The major distinction between Kramers’ theory and the transition state theory lies in the fact that the pre-exponential factor of the latter does not contain a diffusion coefficient and is simply $2\pi /\omega$, where $\omega$ can represent either ${\omega }_o$ or ${\omega }_c$.

According to Equations (3), (4), (5) and (6), the ratios of $t_{cl}/t_{tp}^{c\to o}$ and $t_{op}/t_{tp}^{o\to c}$ are dependent only on ${\omega }^{*}/{\omega }_c$, ${\omega }^{*}/{\omega }_o$ and the height of the free-energy barrier, but are independent of the diffusion coefficient $D^{*}$. As shown in Figure 5b, $t_{op}/t_{tp}^{o\to c}$ is nearly temperature-independent. Given that ${\omega }^{*}/{\omega }_o$ is constant, similar to protein folding [52], the free-energy barrier for the transition from the open state to the closed state should be temperature-dependent, implying that $∆G_o\propto T$. By fitting the two curves, the free-energy barrier of the base pair $∆G_c$ from the closed state to the open state is calculated to be 6.83 kcal/mol, which aligns with the magnitude of enthalpy change $∆H$ from the open state to the closed state. Based on the free-energy change between the open and the closed states, $∆\mathrm{G}=∆G_c-∆G_o=∆H-T∆S$, the free-energy barrier for opening the base pair from the closed state is $∆G_o=T∆S$, where $∆S$ is the entropy change between the closed and open states. This indicates that breaking a single base pair is primarily driven by an increase in enthalpy, ΔH, due to the disruption of base-stacking interactions and hydrogen bonding. Conversely, forming a single base pair mainly involves the entropy loss, ΔS, arising from factors such as the conformation restriction required for pairing with each other, solvent viscosity, and other influences.

3. Materials and Methods

The initial structure for the 6 bp helix with the sequence $\left(\begin{array}{c}{5}^{\prime }-\mathrm{d}\mathrm{A}\mathrm{A}\mathrm{G}\mathrm{A}\mathrm{G}\mathrm{A}-3^{\prime }\\ 3^{\prime }-\mathrm{r}\mathrm{U}\mathrm{U}\mathrm{C}\mathrm{U}\mathrm{C}\mathrm{U}-5^{\prime }\end{array}\right)$ was obtained from its crystal structure (PDB ID: 1DRR). The structure was immersed in a TIP3P [57,58] water box with a water shell of 15 Å in a triclinic box measuring 5.5 × 5.3 × 5.3 nm. The system was subsequently neutralized by the addition of 12 Na⁺ counterions, and additional Na⁺ and Cl⁻ ions were randomly placed to maintain 0.5 M NaCl concentration. The final system included approximately 15,326 atoms. The solvent and counterions were relaxed by energy minimization to equilibrium for 60 ns at 298 K. One of the equilibrium structures was then selected as the starting structure for further simulations at high temperatures.

All simulations were performed using the GROMACS 4.6.7 simulation package [59]. The amber force fields for DNA and RNA were Amber bsc1 and χOL3 [60,61], respectively. Periodic boundary conditions were employed. Temperature coupling was handled by velocity rescaling [62], and pressure coupling used the Parrinello–Rahman barostat algorithm [63]. Electrostatic interactions were evaluated via the Particle-Mesh Ewald (PME) method [64] with a 10 Å direct space cutoff, while Lennard-Jones interactions [65] were truncated at 10 Å. Bond lengths involving hydrogen atoms were constrained using the LINCS algorithm [66], and water molecules were kept rigid using the SETTLE algorithm [67]. The neighboring grid search method was applied and updated every 10 steps. The equations of motion were integrated using the Verlet algorithm with a 2 fs time step, and the coordinates were saved every 2 ps. To specifically investigate the thermodynamic and kinetic properties of the opening and closing of the terminal dA-rU base pair with the nearest neighbors and next-nearest neighbors $\left(\begin{array}{c}{5}^{\prime }-\mathrm{d}\mathrm{A}\mathrm{A}\mathrm{G}\cdots \\ 3^{\prime }-\mathrm{r}\mathrm{U}\mathrm{U}\mathrm{C}\cdots \end{array}\right)$ in the DNA-RNA hybrid, all atoms of DNA-RNA hybrid duplexes except the two terminal nucleotides 5′-dA and 3′-rU were fixed using a harmonic potential with a force constant of 2000 kJ/mol·nm⁻². To accelerate equilibration, the simulation temperatures were set at 370 K, 380 K, 390 K, and 400 K, where the population of the open and closed states each occupied approximately 50%.

4. Conclusions

In summary, through microsecond (4 μs) all-atom molecular dynamics (MD) simulations conducted at various temperatures, the thermodynamic and kinetic parameters of the dA-rU base pair in the DNA-RNA hybrid duplex were derived. Our calculations showed that the entropy change ∆S and the enthalpy change ∆H for opening/closing the dA-rU base pair were consistent with the nearest-neighbor model. These results were close to the thermodynamic parameters of the AU base pair in RNA, aligning well with the experimentally measured single base pair stabilities in RNA and DNA-RNA hybrid duplexes. The average lifetimes of the open state and the closed state showed different dependence on temperature, while the transition path times exhibited only weak temperature dependence. The breaking of a single base pair results from the enthalpy increase ∆H caused by the disruption of the base-stacking interactions and hydrogen bonds. Conversely, the formation of a single base pair is driven primarily by the unfavorable entropy loss ∆S due to accompanying restrictions, such as the torsional angles and solvent viscosity. In conclusion, our simulations provide an atomic-level insight into the opening–closing transition of the dA-rU base pair within a DNA-RNA hybrid duplex and reveal the microscopic mechanism governing the behavior of individual dA-rU base pairs.