Effect of temperature on anisotropic bending elasticity of dsRNA: an all-atom molecular dynamics simulation

Abstract

Employing all-atom molecular dynamics simulations, we examined the temperature-dependent behavior of bending elasticity in double-stranded RNA (dsRNA). Specifically, we focused on the bending persistence length and its constituent components, namely, the tilt and roll stiffness. Our results revealed a near-linear decrease in these stiffness components as a function of temperature, thereby highlighting the increased flexibility of dsRNA at elevated temperatures. Furthermore, our data revealed a significant anisotropy in dsRNA bending elasticity, which diminished with increasing temperature, attributable to marked disparities in tilt and roll stiffness components. We delineated the underlying biophysical mechanisms and corroborated our findings with extant literature. These observations offer salient implications for advancing our understanding of nucleic acid elasticity, and are pertinent to potential medical applications.

---

1 Introduction

RNAs are pivotal in various biological processes, including but not limited to the encapsulation of double-stranded RNA (dsRNA) viral genomes into capsids and structural rearrangements of ribosomes during translation.^1–3 Within these biological contexts, RNA undergoes diverse mechanical deformations in response to alterations in physiological conditions such as stretching, bending, and twisting. Concurrently, dsRNA has garnered attention as a promising nanomaterial for biological and nanomedical applications.³ Bend elasticity is a fundamental mechanical property of nucleic acids and is instrumental in shaping their three-dimensional structures under physiological conditions.⁴ Accordingly, elucidating the impact of environmental variables such as ionic strength and temperature on dsRNA bend elasticity is crucial.

Traditionally, bending elasticity is quantified either via the bending persistence length or bending stiffness. For dsDNA, the bending persistence length is approximately 50 nm under standard physiological conditions, as described by the wormlike chain (WLC) model.^5,6 Notably, the bending elasticity of dsDNA is sensitive to external conditions, including ionic concentration^7–11 and temperature.^12–21 Temperature elevations typically manifest as reductions in the bending persistence length of dsDNA, which is attributable to thermal fluctuations within the molecular milieu.²² For examples, the experiments suggested that the bending persistence length decreased as the temperature was increased,¹⁷ which have been reproduced in Monte Carlo (MC) and all-atom molecular dynamics (MD) simulations based on various models.^19,20,23

For dsRNA, both experimental and computational data suggest a bending persistence length of approximately 60 nm under normal physiological conditions, surpassing that of dsDNA.^24–26 For examples, the magnetic tweezer (MT) experiments reported that the dsRNA possesses the bending persistence lengths of l_B = 60 ± 1 nm,²⁴ l_B = 57 ± 2 nm (ref. 25) and l_B = ∼61 nm.²⁶ The all-atom MD data reproduced that the dsRNA have the bending persistence length of l_B = 69 ± 4 nm,²⁷ l_B = 66.3 nm (ref. 25) and l_B = 66.99 ± 1.38 nm.²⁸ The observed disparities between the MT and MD data can be attributed to variations in experimental conditions, such as ionic strength and sequence composition. The influence of sequence composition on elasticity mainly refers to the length and order of base pairs. Generally, the sequence length used in experiments is several thousand bp, while the sequence length used in all-atom simulations is several tens of bp due to the limitations of computational performance. It is noteworthy that the bending persistence length of dsRNA diminishes with an increase in ion concentration.²⁶ Nevertheless, the impact of temperature on the bending persistence length of dsRNA remains an open avenue for investigation.

In the extant Marko and Siggia (MS) elastic model, bending elasticity is partitioned into two orthogonal modes, encapsulated by the tilt and roll angles.²⁹ This partitioning underpins the observed anisotropy in bending rigidities, manifesting as differential flexibility along the major and minor grooves of nucleic acid helices.³⁰ The MS model, extensively employed in dsDNA studies, accounts for the observed bending anisotropies and associated twist–bend couplings.^31–35 Previous research cites tilt stiffness values of approximately 100 nm and roll stiffness values of approximately 39 nm for dsDNA at room temperature.^31,32,36,37 The inherent asymmetry between these values correlates with the observed molecular geometry of the DNA.³⁶ Such anisotropic behavior has also been observed in dsRNA, albeit with sequence-dependent variations in the extent of anisotropy.³⁸ Recently, Dohnalová et al. used single-molecule MT measurements to determine the temperature-dependence of the dsRNA twist and observed that dsRNA unwinds with increasing temperature, which was correctly predicted by the all-atom MD simulations.³⁹ Dabin et al. investigated the thermal denaturation of dsRNA using the atomistic simulations by varying the temperature in a wide range, in which the sequences and force fields were considered.⁴⁰

In the present study, we employed all-atom molecular dynamics simulations of a representative dsRNA sequence to scrutinize the temperature-dependent properties of bending elasticity. Our study concentrates on the temperature effects on the bending persistence length, tilt, and roll stiffness, as well as the bending anisotropy, guided by the WLC and MS models. Section 2 outlines the simulation methodology, Section 3 discusses temperature-dependent findings, and Section 4 concludes.

2 Simulation models and method

2.1 All-atom MD simulation

In the current simulations, we selected the initial A-type dsRNA with a 16 bp sequence of 5′-GCGC AAUG GAGU ACGC-3′, which has been used in the previous works.^41,42 Our all-atom MD simulations are similar to those in previous studies,^21,43,44 and we only briefly describe the content involved in the current simulations. We constructed the initial structure file of the RNA using the UCSF Chimera 1.15 software application,⁴⁵ as shown in Fig. 1(a), where there are approximately 1026 atoms in the 16 bps sequence. The dsRNA sequence was immersed in an 8 × 8 × 8 nm³ simulation box under periodic boundary conditions. Subsequently, NaCl⁴⁶ and water molecules were added corresponding to 100 mM salt concentration, where the TIP3P model⁴⁷ was used to define the water structure. We added Na⁺ and Cl⁻ ions using by Joung and Cheatham's ion model⁴⁸ to the environment to maintain a constant concentration. Because nucleic acids are negatively charged, additional 30 Na⁺ ions are dumped into the solution, which neutralizes the negative charge in the solution environment.

Fig. 1 (a) Representative diagrams for the molecular structure of 16 bps dsRNA with 5′-GCGC AAUG GAGU ACGC-3′ sequence. (b) Root mean square deviation (RMSD) curves of the 10 bases fragment in the center of the dsRNA at T = 300 K, where the black line indicates the average value of the relevant parameter every 2 ns, which reached equilibrium after 100 ns. (c) Sketch diagrams for the tilt angle, roll angle, and twist angle on the molecular structure.

Our MD simulations were performed using the Gromacs 4.6 software with the Amber OL3 force fields.^49–53 Before the MD simulations, the system needs to be pre-equilibrated to ensure that the MD simulation is performed in an isothermal isobaric environment.^54,55 We started with an energy minimization system with a restraining force and subsequently an NVT and NPT simulations. In NVT, the system is warmed up to the desired temperature T by the V-rescale method.^54,55 Specifically, this scheme uses the velocity-rescaling method with random terms for temperature coupling to produce the NVT canonical ensemble, retaining the Berendsen thermostat's advantages, which have the first-order decay of the temperature deviation and the absence of oscillations.^54,55 After the system was preheated to the desired temperature T, the system pressure was adjusted to 1 atm using the NPT simulation. When these pre-equilibration processes were completed, and a 1000 ns MD simulation was performed at a fixed temperature T without the restraining force. We adjusted the temperature T from 280 K to 320 K, with a step of ΔT = 10 K, where the temperature regulation was similar to those of all-atom MD simulations with Amber force fields.^{16,20,21,56,57} We checked the validity of the simulation by calculating the root mean square deviation (RMSD) of the system. An example is shown in the RMSD plots for 16 bps dsRNA in Fig. 1(b). The RMSD indicates the extent of structural changes in the dsRNA molecule as follows:

where δ_i is the displacement of the i atom from moment 0 to moment t with a time step of Δt = 2 ps. We take each 2 ns as the average value and draw the black line in the middle as shown in Fig. 1(b). We also provided more RMSD plots at various temperatures in Fig. S1 of ESI.† These RMSD data indicated that these systems with various temperatures reach to the equilibrium states after 100 ns, respectively. In all simulations, elasticity data from 100 ns to 1000 ns were used for the subsequent statistical analysis. To avoid end effects, we removed 3 bps from each end of the sequence in all data analysis and only selected the central 10 bps to analyze the overall elastic properties of dsRNA.

2.2 WLC and MS models

In this study, we focused on the elastic parameters of the MS model and WLC models. In WLC model, the bend elasticity can be described as follows^9,58–60Here, l_B is the bending persistence length, p(θ, l) is the probability distribution; θ is the bending angle formed by a dsRNA spanning 6 bps;⁶¹ l is a constant representing the average 6 bps profile length. The l_B value was obtained by fitting the quadratic function to eqn (2) using the MD simulation data. Here, we note that by removing several bps from each helix end in the data analysis, the short nucleic acid sequence can also be described by WLC model with the bending persistence length using eqn (2).⁶²

The MS model was originally proposed by Marko and Siggia,²⁹ where three angles Ω₁, Ω₂ and Ω₃ was used to describe the bend and twist elasticities of DNA in recent years.^31,34,35,63 These three rotational angles, Ω₁, Ω₂ and Ω₃, are related to the tilt, roll and twist in the rigid base-pair,³⁷ as shown in Fig. 1(c). In the MS model, by considering the symmetry and lowest order of Ω, one can express the energy functional as follows²⁹

where k_B is the Boltzmann constant, T is the temperature and the dots denote the higher-order terms. The A₁ and A₂ are two bending stiffnesses, C is the torsional stiffness and G is the twist–bend coupling constant. These elastic parameters have length dimensions; thus A₁, A₂ and C also denote the corresponding persistence lengths. The MS model can be reduced to the twistable worm-like chain (TWLC) model by neglecting the twist–bend coupling G, and taking the energy functional as follows³²

One can further reduce eqn (4) into WLC model by ignoring Ω₃.

In the MS model, A₁, A₂, C and G are related to the elastic moduli K_ττ, K_ρρ, K_ωω and K_ρω multiplied by L₀ in unit of k_BT, where L₀ is the contour length of dsRNA. According to the random thermal fluctuation, the elastic modulus K can be expressed as²²

Here, matrix C is a covariance matrix with the elements

cov(i,j) = c_ijσ_iσ_j (6)

where c_ij is the Pearson coefficient of i and j, σ_i and σ_j are the standard deviations of i and j. The i and j denote the cumulative tilt τ, roll ρ and twist ω, which are MD-associated parameters in the simulations.

Generally, three angular parameters, tilt, roll and twist, can be used to describe the bending and twisting elasticities of dsRNA, as well as its coupling, in the rigid base-pair (RBP) model, as shown in eqn (5). Here, we try to describe the relationship between MS, TWLC and WLC models, starting from the RBP model. In MS model, only the coupling between the roll and twist, is considered, since the tilt is weakly coupled to roll and twist.²⁹ However, the TWLC model further ignore the coupling between roll and twist angles.³² In WLC model, it only considers the tilt and roll angles by ignoring the twist angle, which can be used to describe the bending elasticities. Such approximation of the WLC model makes the MS model closer to the experiments.

3 Results and discussion

In this study, we focus on the influence of temperature on the bend elasticity of dsRNA, which is characterized by the bending persistence length based on the WLC model and the bending anisotropy based on the MS model. The temperature was adjusted from T = 280 K to T = 320 K with steps of ΔT = 10 K to demonstrate the temperature-dependent bend elasticity of dsRNA. In Subsection 3.1, we discuss the dependence of bending elasticities in Fig. 2 and 3; we analyzed the temperature dependence of tilt and roll stiffnesses, as well as their bending anisotropy, based on the MS model in Fig. 4–7 in Subsection 3.2.

3.1 Temperature-dependent bend elasticity

We plotted the bending persistence length l_B as a function of temperature T in Fig. 2, according to eqn (1) based on the WLC model. In Fig. 2(a), a typical example is shown, where the is plotted as a function of the bending angle θ at T = 300 K. We used the quadratic curve to fit the MD data and obtained the bending persistence length l_B = 62.06 ± 0.41 nm at T = 300 K. Our MD data are in line with the available all-atom MD simulation data where l_B ranges from 66 nm to 69 nm.^25,27,28 Previous experiments have suggested that the dsRNA has a persistence length of approximately l_B = 57–64 nm under the various salt concentrations.^24–26 In particular, Abels et al. used two single-molecule techniques, the MT and AFM, to measure the bending persistence length of dsRNA at room temperature (T = ∼298 K) and obtained l_B = 63.8 ± 0.7 nm by MT and l_B = 62 ± 2 nm by AFM.⁶⁴ We listed these data in detail in Table S3 of ESI† where the sequence length and ion concentration were also listed. Actually, the bending persistence length of dsRNA decreases as the concentration of CoHex³⁺ increases¹¹ and also has a sequence-dependence.³⁸ Thus, the main deviation between the current MD results and experimental data was probably due to the different salt concentrations and sequence lengths used in the current simulation and previous works. For a convenient comparison, we inserted these available data into Fig. 2(b), where the bending persistence lengths were plotted as a function of temperature. Here, we predicted the temperature dependence of bending persistence length l_B, as shown in Fig. 2(b), where the temperature T over a range from T = 280 K to T = 320 K. The bending persistence lengths at other temperatures are also obtained by the same fitting methods as those at T = 300 K [see Fig. S2 of ESI†]. We observed a linear relationship between the bending persistence length l_B and temperature T, with the slope of k_l = −0.342 nm K⁻¹. Previous experiments have reported the linear dependence of temperature on the bending persistence lengths of dsDNA,^14,17 and all-atom simulations have also reproduced the linear temperature dependence on the bending persistence lengths of dsDNA where the bending persistence length l_B decreases as temperature T increases, with a slope of −0.29 nm K⁻¹.²¹

Fig. 2 (a) The relationship between and bending angle θ at T = 300 K. The bending angle θ is formed by six consecutive base pairs on each of the 10 base segments at the center of the dsRNA. (b) Temperature dependence of bending persistence length of dsRNA. The line is a fitting result with a slope of k_l = −0.342 nm K⁻¹. The available data were also inserted for convenient comparison.

We used the rotational parameters to understand the temperature dependence of l_B, as shown in Fig. 3. We plotted the bending angle 〈θ〉 as function of temperature T in Fig. 3(a), where θ is bending angle between two successive bps and 〈⋯〉 denotes assemble average. We obtained the 〈θ〉 from the MD trajectory data, and 〈θ〉 = 9.05° at T = 300 K. We note that the average bending angle 〈θ〉 is much greater than that of dsDNA.^21,27,62 Furthermore, we predicted the linear dependence of the average bending angle 〈θ〉 on the temperature T, as shown in Fig. 3(a). The results showed that the bending angle 〈θ〉 increased as the temperature T increased, with a slope of k_θ = 0.021° K⁻¹, which is similar to the dsDNA case.²¹ We are aware that the smaller bending angle reflects the stronger dsDNA rigidity, in which the upward trend illustrated that the dsRNA chain becomes more flexible with the increasing T, supported by the short WLC model about the conformational variations.^58,65 In the TWLC model, the bending angle is related to two orthogonal angles, the tilt and roll angles, which hold θ² = τ² + ρ².^31,35 To analyze the bending elasticity in more detail, we plotted the 〈τ〉 and 〈ρ〉 as functions of temperature T, as shown in Fig. 3(b) and (c), where τ and ρ are the tilt and roll angles, respectively, between two successive bps. We also listed the more detailed data about the tilt and roll angles at various temperatures in Table S1 of ESI.† We obtained 〈τ〉 = 0.08° and 〈ρ〉 = 9.05° at T = 300 K, which are in line with the available data from the various dsRNA sequences,^25,27,28 for example, 〈ρ〉 = 7.5 ± 3.47° and 〈τ〉 = −0.01 ± 4.36° for a 16 bps dsRNA sequence used in the previous MD simulation.²⁸ The results showed that the average roll angle 〈ρ〉 has a linear relationship with temperature T on a slope of k_ρ = 0.021° K⁻¹. However, the average tilt angle 〈τ〉almost maintain unchanged value, with a slope of k_τ = 0.000° K⁻¹, as shown in Fig. 3(b) and (c), which suggests that the bending elasticity is mainly derived from the roll angles.⁶⁶

Fig. 3 (a) The temperature dependent bending angle θ. The line is a fitting result with a slope of k_θ = 0.021° K⁻¹. (b) The temperature dependent roll angle ρ. The lines are fitting results with slope of k_ρ = 0.021° K⁻¹. (c) The temperature dependent tilt angle τ. The lines are fitting results with slope of k_τ = 0.000° K⁻¹.

3.2 Temperature-dependent tilt and roll stiffnesses

Tilt, roll, and twist are three rotational parameters that play important roles in the bend and twist stiffnesses. In the MS model, the tilt and roll are the two orthogonal components of the bend elasticity. In the current simulations, we concentrated on the bend elasticities involving the tilt and roll angles, as well as the bending anisotropy in Fig. 4 and 5.

Fig. 4 (a) An example of correlation between tilt (τ) and roll (ρ). The data correction coefficient is c_τρ = −0.029. (b) Temperature dependence of correlation coefficient c_τρ, c_τω, and c_ρω on temperature T. The line is a fitting result with a slope of c_kτρ = 2.900 × 10⁻⁴ K⁻¹, c_kτω = −2.000 × 10⁻⁵ K⁻¹, and c_kρω = −8.500 × 10⁻⁴ K⁻¹.

Fig. 5 The temperature dependence of tilt stiffness A₁, roll stiffness A₂. (a) The function of tilt stiffness A₁ as a function of temperature T, and the line is a fitting result with a slope of k₁ = −0.580 nm K⁻¹. (b) The function of roll stiffness A₂ as a function of temperature T, and the line is a fitting result with a slope of k₂ = −0.167 nm K⁻¹.

First, we analyzed the tilt and roll correlations and showed an example at T = 300 K in Fig. 4(a), in which the data between the tilt and roll have negative correlations, with a negative slope of c_τρ = −0.029. This weak correlation indicates that the tilt angle was almost independent of the roll angle at T = 300 K. A previous MD simulation also indicated that there was a weak correlation between the tilt and roll angles for dsRNA and dsDNA.²⁵ Correlations between tilt and roll are presented in a manner similar to the previous simulations.^43,57 We also provided the tilt and roll correlations at other temperature in Fig. S3 of ESI,† which indicated the correlations between the tilt and roll are also weak. To demonstrate the temperature-dependence of this correlation clearly, we plotted the correlation coefficient as a function of the temperature T, as shown in Fig. 4(b). We then observed a linear relationship between the correlation coefficient c_τρ and temperature T with a slope of c_kτρ = 2.900 × 10⁻⁴ K⁻¹. The very weak upward trend indicates that the correlation between the tilt and roll angles is independent of the temperature T. Actually, the correlations of tilt and roll to the third angle, twist, cannot be neglected in dsRNA, as shown in Fig. 4(b) [see the tilt–twist correlations and roll–twist correlations at various temperatures in Fig. S4 and S5 of ESI†], which differs from the available data for dsDNA.²⁵ This non-correlation between the tilt and roll enables us to analyze the bend elasticity based on the tilt and roll stiffnesses.

We then considered two elastic components, tilt and roll stiffnesses, for the bend elasticity illustrated in the MS model. The MS model has been successfully used in dsDNA to explain the bend and twist elasticities by considering the twist-bending coupling.^31–35 This is due to undistortion of dsDNA, where the bending angle is relatively small, which leads to the rotational symmetry of tilt angle. The dsRNA is more bent than dsDNA where the average roll between two successive bps is 〈ρ〉 = 9.05°. This causes the nonzero terms K_τω ≠ 0 [refer to the K elements at various temperatures in Table S2 of ESI†], leading to the deviation of the MS model. However, the K_ρω is much larger than K_τω, which enables us to analyze the tilt and roll components for bending elasticities. This is also supported by the fact that the Pearson correlation coefficients between the roll and twist are much larger than those between tilt and twist [refer to Table S2 in ESI† for K_ρω and K_τω at various temperatures]. We obtained the tilt and roll modulus, K_ττ and K_ρρ, according to eqn (5) and (6), and plotted the tilt and roll stiffnesses (tilt and roll persistence lengths), A₁ and A₂, as function of temperature T, in Fig. 5. Here, the tilt stiffness A₁ is determined by^31,32

where L₀ denotes the contour length. We obtained A₁ = 120.80 nm by substituting K_ττ = 169.39 pN nm and A₂ = 22.46 nm by substituting K_ρρ = 31.50 pN nm at T = 300 K, where L₀ = 2.954 nm. The available data showed that the tilt stiffness A₁ are about ∼100 nm for dsDNA,^31,32,36,37 whereas the roll stiffness A₂ = 39 nm (ref. 37) and A₂ = 38.8 nm(ref. 32) for dsDNA near the room temperature. Our observations showed that dsRNA has a smaller roll stiffness A₂ than those in dsDNA. It is meaningful to use the asymmetric parameter to describe the bending anisotropy, which is defined as follows^29,30

B = (A₁ − A₂)/2 (8)

In the current simulations, the dsRNA has an asymmetry parameter of B = 49.17 nm, which indicates a much larger bending anisotropic bending elasticity. This asymmetry parameter was much larger than that in dsDNA where the asymmetry parameter was suggested to be B = 19 nm.³⁰ This anisotropic bending elasticity was also observed in the MD simulations of the dsRNA and dsDNA.³⁸ Importantly, our simulation results predicted that the tilt stiffness A₁ and roll stiffness A₂ are temperature-dependent, as shown in Fig. 5(a) and (b), respectively. The data fitting indicated a linear relationship between A₁ and T with a slope of k₁ = −0.580 nm K⁻¹, while the same is true for the roll angle where A₂ decreases with T by a slope of k₂ = −0.167 nm K⁻¹. Our observations were supported by a theoretical analysis in which there was a linear relationship between the stiffness matrix K and temperature T.¹⁶ We also show the effect of temperature on the anisotropic bending elasticity in Fig. S6.† Our simulation results indicate that asymmetry parameter B decreases from 53.85 nm to 44.42 nm when the temperature T increases from 280 K to 320 K, with a slope of k₃ = −0.208 nm K⁻¹, indicating that the bending anisotropy become weaker as temperature increases.

3.3 Thermal fluctuations in tilt and roll angles

According to thermal fluctuations, the elasticity can be described by the covariances between the deformation variables.²² Then, we delve into the mechanism governing the temperature dependence of tilt stiffness A₁ and roll stiffness A₂, as shown in Fig. 6. eqn (5) and (6) suggest that the stiffness matrix K can be described by the Pearson correlation coefficient cov(i, j) and thermal fluctuation σ_i². This enables us to estimate that the tilt stiffness A₁ and roll stiffness A₂ have approximate forms, similar to the stretch–twist elastic matrix,^21,22

Fig. 6 (a) The variance of cumulative tilt σ_τ² as a function of temperature T. The line is a fitting result with a slope of k_στ = 1.080 × 10⁻⁴ rad² K⁻¹. (b) The variance of cumulative roll σ_ρ² as a function of temperature T. The line is a fitting result with a slope of k_σρ = 9.310 × 10⁻⁴ rad² K⁻¹.

As the contour length L₀ remains unchanged values with temperature T, we present the factors σ_τ² and σ_ρ² as functions of temperature T to elucidate the decline in tilt stiffness A₁ and roll stiffness A₂, as shown in Fig. 6. The results indicate that both σ_τ² and σ_ρ² increase as the temperature T increases, displaying nearly linear trends, which is similar to the data obtained from MD simulations of dsDNA.²¹ These findings suggest that the reductions in tilt stiffness A₁ and roll stiffness A₂ are attributed to the thermal fluctuation components σ_τ² and σ_ρ², respectively. However, the fitting results suggest that the decline rates differ between the σ_τ² and σ_ρ² as the temperature T increases. Specifically, σ_τ² increases with a slope of k_στ = 1.080 × 10⁻⁴ rad² K⁻¹ while k_σρ = 9.310 × 10⁻⁴ rad² K⁻¹ for σ_ρ². This implies that the thermal fluctuation of the roll angle σ_ρ² depends more obviously on the temperature T than on the tilt angle, which plays a stronger role in the T-dependent roll stiffness A₂.

To understand the thermal fluctuations in more detail, we investigated the thermal fluctuation components σ_t² and σ_r² for each base pair. An example at T = 300 K was shown in Figs. 7(a) and (b), respectively. The data showed that σ_r² is greater than σ_t² in each bp, and generally the σ_t² and σ_r² in U and A bps are greater than those in C and G bps, which is in line with the sequence-dependent data.^38,67 The min values of σ_t² and σ_r² appear in the center of the sequence and the larger values at the two ends, exhibiting the flexibilities of these two ends. These flexibilities are consistent with the flexibility origin of short sequences.^62,68 Here, we used the sequence-dependent σ_t² and σ_r² to clearly illustrate the anisotropic bending elasticities originating from the thermal fluctuations.

Fig. 7 (a) The variance of tilt σ_t² for each bp T = 300 K. (b) The variance of roll σ_r² for each bp T = 300 K. The dashed curves denote the varying trends along dsRNA sequences.

4 Conclusions

In the present study, we used all-atom MD simulations to investigate the effects of temperature on the bending elasticity of dsRNA. We chose a short dsRNA sequence of 16 bps to analyze the bending persistence length, tilt, and roll stiffness of the dsRNA. We concentrated on the influences of temperature on bending elasticities and bending anisotropy based on WLC and MS models, and analyzed the mechanism of temperature dependence by thermal fluctuations from the rotational deformation variables.

In the current simulations, we obtained the bending persistence length l_B = 62.06 ± 0.41 nm for dsRNA at T = 300 K, which is in line with the available experimental and MD simulation data. We predicted a linear relationship between the bending persistence length l_B and temperature T, with a slope of k_l = −0.342 nm K⁻¹. This linear relationship was similar to that observed for dsDNA. To analyze the bending elasticity in more detail, we investigated the 〈τ〉 and 〈ρ〉 as functions of temperature T, and the results showed that the average roll angle 〈ρ〉 has a linear relationship with temperature T at a slope of k_ρ = 0.021° K⁻¹, while the average tilt angle 〈τ〉 remained almost unchanged, with a slope of k_τ = 0.000° K⁻¹, which suggests that the bending elasticity is mainly from the roll angles. We then investigated the tilt and roll stiffness as well as its bending anisotropy. We obtained a tilt stiffness A₁ = 120.80 nm and roll stiffness A₂ = 22.46 nm at T = 300 K. DsRNA has an asymmetry parameter of B = 49.17 nm, which indicates a much larger bending anisotropic bending elasticity than that of dsDNA, where the asymmetry parameter was suggested to be B = 19 nm. We predicted that the tilt stiffness A₁ and roll stiffness A₂ decrease linearly with increasing temperature T in linear manners, with a slope of k₁ = −0.580 nm K⁻¹ for tilt stiffness, while the same is true for the roll angle, where A₂ decreases with T by a slope of k₂ = −0.167 nm K⁻¹. Our MD simulation results showed that asymmetry parameter B decreases from 53.85 nm to 44.42 nm when the temperature T increases from 280 K to 320 K, with a slope of k₃ = −0.208 nm K⁻¹. This suggests that the bending anisotropy weakened as the temperature increased. Then, we delve into the mechanism governing the temperature dependence of tilt stiffness A₁ and roll stiffness A₂, which suggests that the thermal fluctuation of the roll angle σ_ρ² depends more on the temperature T than on the tilt angle. As we know, the bending of RNA plays a key role in the transmission of genetic information and organisms survive, as well as the nucleic acid nanostructures involved in functions related to drug delivery. Our all-atom MD simulations provide a deeper understanding of the bending elasticity of dsRNA, which probably has potential applications in these aspects.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank for the financial supports from the Program of National Natural Science Foundation of China (Grant No. 12374215 and 11875205).

Notes and references

1. F. Tama, M. Valle, J. Frank and C. L. Brooks III, Proc. Natl. Acad. Sci. U. S. A., 2003, 100, 9319–9323.
2. B. S. Schuwirth, M. A. Borovinskaya, C. W. Hau, W. Zhang, A. Vila-Sanjurjo, J. M. Holton and J. H. D. Cate, Science, 2005, 310, 827–834.
3. Q. Vicens and J. S. Kieft, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2112677119.
4. A. Basu, D. G. Bobrovnikov and T. Ha, J. Mol. Biol., 2021, 433, 166861.
5. A. Aggarwal, S. Naskar, A. K. Sahoo, S. Mogurampelly, A. Garai and P. K. Maiti, Curr. Opin. Struct. Biol., 2020, 64, 42–50.
6. A. Marin-Gonzalez, J. Vilhena, R. Perez and F. Moreno-Herrero, Q. Rev. Biophys., 2021, 54, e8.
7. A. Garai, S. Saurabh, Y. Lansac and P. K. Maiti, J. Phys. Chem. B, 2015, 119, 11146–11156.
8. A. Brunet, C. Tardin, L. Salomé, P. Rousseau, N. Destainville and M. Manghi, Macromolecules, 2015, 48, 3641–3652.
9. A. V. Drozdetski, I. S. Tolokh, L. Pollack, N. Baker and A. V. Onufriev, Phys. Rev. Lett., 2016, 117, 028101.
10. S. Guilbaud, L. Salomé, N. Destainville, M. Manghi and C. Tardin, Phys. Rev. Lett., 2019, 122, 028102.
11. H. Fu, C. Zhang, X.-W. Qiang, Y.-J. Yang, L. Dai, Z.-J. Tan and X.-H. Zhang, Phys. Rev. Lett., 2020, 124, 058101.
12. J.-S. Park, K. J. Lee, S.-C. Hong and J.-Y. Hyon, J. Korean Phys. Soc., 2008, 52, 1927–1931.
13. R. A. Forties, R. Bundschuh and M. G. Poirier, Nucleic Acids Res., 2009, 37, 4580–4586.
14. S. Geggier, A. Kotlyar and A. Vologodskii, Nucleic Acids Res., 2011, 39, 1419–1426.
15. N. Theodorakopoulos and M. Peyrard, Phys. Rev. Lett., 2012, 108, 078104.
16. S. Meyer, D. Jost, N. Theodorakopoulos, M. Peyrard, R. Lavery and R. Everaers, Biophys. J., 2013, 105, 1904–1914.
17. R. P. Driessen, G. Sitters, N. Laurens, G. F. Moolenaar, G. J. Wuite, N. Goosen and R. T. Dame, Biochemistry, 2014, 53, 6430–6438.
18. A. Brunet, L. Salomé, P. Rousseau, N. Destainville, M. Manghi and C. Tardin, Nucleic Acids Res., 2018, 46, 2074–2081.
19. J. M. Schurr, Biochemistry, 2019, 251, 106146.
20. H. Dohnalová, T. Drsata, J. Sponer, M. Zacharias, J. Lipfert and F. Lankas, J. Chem. Theory Comput., 2020, 16, 2857–2863.
21. Y. Zhang, L. He and S. Li, J. Chem. Phys., 2023, 158, 094902.
22. F. Lankaš, O. Gonzalez, L. Heffler, G. Stoll, M. Moakher and J. H. Maddocks, Phys. Chem. Chem. Phys., 2009, 11, 10565–10588.
23. J. J. Delrow, P. J. Heath and J. M. Schurr, Biophys. J., 1997, 73, 2688–2701.
24. E. Herrero-Galan, M. E. Fuentes-Perez, C. Carrasco, J. M. Valpuesta, J. L. Carrascosa, F. Moreno-Herrero and J. R. Arias-Gonzalez, J. Am. Chem. Soc., 2013, 135, 122–131.
25. J. Lipfert, G. M. Skinner, J. M. Keegstra, T. Hensgens, T. Jager, D. Dulin, M. Köber, Z. Yu, S. P. Donkers and F.-C. Chou, et al., Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 15408–15413.
26. C. Zhang, H. Fu, Y. Yang, E. Zhou, Z. Tan, H. You and X. Zhang, Biophys. J., 2019, 116, 196–204.
27. Y. Zhang, M. Yan, T. Huang and X. Wang, Adv. Theory Simul., 2023, 6, 2200534.
28. K. B. Chhetri, A. Sharma, S. Naskar and P. K. Maiti, Nanoscale, 2022, 14, 6620–6635.
29. J. F. Marko and E. D. Siggia, Macromolecules, 1994, 27, 981–988.
30. D. Norouzi, F. Mohammad-Rafiee and R. Golestanian, Phys. Rev. Lett., 2008, 101, 168103.
31. S. K. Nomidis, F. Kriegel, W. Vanderlinden, J. Lipfert and E. Carlon, Phys. Rev. Lett., 2017, 118, 217801.
32. E. Skoruppa, M. Laleman, S. K. Nomidis and E. Carlon, J. Chem. Phys., 2017, 146, 214902.
33. E. Skoruppa, S. K. Nomidis, J. F. Marko and E. Carlon, Phys. Rev. Lett., 2018, 121, 088101.
34. M. Caraglio, E. Skoruppa and E. Carlon, J. Chem. Phys., 2019, 150, 135101.
35. S. K. Nomidis, E. Skoruppa, E. Carlon and J. F. Marko, Phys. Rev. E, 2019, 99, 032414.
36. F. Mohammad-Rafiee and R. Golestanian, Phys. Rev. Lett., 2005, 94, 238102.
37. Y. A. G. Fosado, F. Landuzzi and T. Sakaue, Phys. Rev. Lett., 2023, 130, 058402.
38. I. Faustino, A. Pérez and M. Orozco, Biophys. J., 2010, 99, 1876–1885.
39. H. Dohnalová, M. Seifert, E. Matouskova, M. Klein, F. S. Papini, J. Lipfert, D. Dulin and F. Lankas, J. Phys. Chem. B, 2024, 128, 664–675.
40. A. Dabin and G. Stirnemann, Biochemistry, 2024, 307, 107167.
41. K. Liebl, T. Drsata, F. Lankas, J. Lipfert and M. Zacharias, Nucleic Acids Res., 2015, 43, 10143–10156.
42. A. Marin-Gonzalez, J. Vilhena, R. Perez and F. Moreno-Herrero, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 7049–7054.
43. X.-W. Qiang, C. Zhang, H.-L. Dong, F.-J. Tian, H. Fu, Y.-J. Yang, L. Dai, X.-H. Zhang and Z.-J. Tan, Phys. Rev. Lett., 2022, 128, 108103.
44. F.-J. Tian, C. Zhang, E. Zhou, H.-L. Dong, Z.-J. Tan, X.-H. Zhang and L. Dai, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2218425120.
45. E. F. Pettersen, T. D. Goddard, C. C. Huang, G. S. Couch, D. M. Greenblatt, E. C. Meng and T. E. Ferrin, J. Comput. Chem., 2004, 25, 1605–1612.
46. L. Bao, X. Zhang, Y.-Z. Shi, Y.-Y. Wu and Z.-J. Tan, Biophys. J., 2017, 112, 1094–1104.
47. S. Izadi, R. Anandakrishnan and A. V. Onufriev, J. Phys. Chem. Lett., 2014, 5, 3863–3871.
48. I. S. Joung and T. E. Cheatham III, J. Phys. Chem. B, 2008, 112, 9020–9041.
49. B. Hess, C. Kutzner, D. Van Der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
50. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
51. R. Galindo-Murillo, J. C. Robertson, M. Zgarbova, J. Sponer, M. Otyepka, P. Jurecka and T. E. Cheatham III, J. Chem. Theory Comput., 2016, 12, 4114–4127.
52. I. Ivani, P. D. Dans, A. Noy, A. Pérez, I. Faustino, A. Hospital, J. Walther, P. Andrio, R. Goñi and A. Balaceanu, et al., Nat. Methods, 2016, 13, 55–58.
53. K. X. Xiong, K. Xi, L. Bao, Z. L. Zhang and Z. J. Tan, Acta Phys. Sin., 2018, 67, 108701.
54. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
55. Q. Ke, X. Gong, S. Liao, C. Duan and L. Li, J. Mol. Liq., 2022, 365, 120116.
56. F. Kriegel, C. Matek, T. Dršata, K. Kulenkampff, S. Tschirpke, M. Zacharias, F. Lankaš and J. Lipfert, Nucleic Acids Res., 2018, 46, 7998–8009.
57. C. Zhang, F. Tian, Y. Lu, B. Yuan, Z.-J. Tan, X.-H. Zhang and L. Dai, Sci. Adv., 2022, 8, eabn1384.
58. A. K. Mazur, Phys. Rev. Lett., 2007, 98, 218102.
59. A. Polley, J. Samuel and S. Sinha, Phys. Rev. E, 2013, 87, 012601.
60. X.-W. Qiang, H.-L. Dong, K.-X. Xiong, W. Zhang and Z.-J. Tan, Commun. Theor. Phys., 2021, 73, 075601.
61. V. Ortiz and J. J. de Pablo, Phys. Rev. Lett., 2011, 106, 238107.
62. Y.-Y. Wu, L. Bao, X. Zhang and Z.-J. Tan, J. Chem. Phys., 2015, 142, 125103.
63. S. Nomidis, M. Caraglio, M. Laleman, K. Phillips, E. Skoruppa and E. Carlon, Phys. Rev. E, 2019, 100, 022402.
64. J. Abels, F. Moreno-Herrero, T. Van der Heijden, C. Dekker and N. H. Dekker, Biophys. J., 2005, 88, 2737–2744.
65. R. Lavery, M. Moakher, J. H. Maddocks, D. Petkeviciute and K. Zakrzewska, Nucleic Acids Res., 2009, 37, 5917–5929.
66. A. Noy and R. Golestanian, Phys. Rev. Lett., 2012, 109, 228101.
67. O. Gonzalez, D. Petkevičiūtė and J. Maddocks, J. Chem. Phys., 2013, 138, 055102.
68. M. Segers, E. Skoruppa, J. A. Stevens, M. Vangilbergen, A. Voorspoels and E. Carlon, J. Chem. Phys., 2021, 155, 027101.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4ra02354d

---