Predictions and analyses of RNA nearest neighbor parameters for modified nucleotides

Abstract

The most popular RNA secondary structure prediction programs utilize free energy (ΔG°37) minimization and rely upon thermodynamic parameters from the nearest neighbor (NN) model. Experimental parameters are derived from a series of optical melting experiments; however, acquiring enough melt data to derive accurate NN parameters with modified base pairs is expensive and time consuming. Given the multitude of known natural modifications and the continuing use and development of unnatural nucleotides, experimentally characterizing all modified NNs is impractical. This dilemma necessitates a computational model that can predict NN thermodynamics where experimental data is scarce or absent. Here, we present a combined molecular dynamics/quantum mechanics protocol that accurately predicts experimental NN ΔG°37 parameters for modified nucleotides with neighboring Watson-Crick base pairs. NN predictions for Watson-Crick and modified base pairs yielded an overall RMSD of 0.32 kcal/mol when compared with experimentally derived parameters. NN predictions involving modified bases without experimental parameters (N6-methyladenosine, 2-aminopurineriboside, and 5-methylcytidine) demonstrated promising agreement with available experimental melt data. This procedure not only yields accurate NN ΔG°37 predictions but also quantifies stacking and hydrogen bonding differences between modified NNs and their canonical counterparts, allowing investigators to identify energetic differences and providing insight into sources of (de)stabilization from nucleotide modifications.

Figure 1.

Base pairs whose nearest neighbor free energy parameters were predicted or utilized for comparison in this work. (A) The Watson–Crick adenosine–uridine base pair (top) and similar modified derivatives: N⁶-methyladenosine·uridine, 2,6-diaminopurineriboside·uridine, and 2-aminopurineriboside·uridine. (B) The Watson-Crick guanosine-cytidine base pair (top) and similar modified derivatives: inosine·cytidine, guanosine·5-methylcytidine, and isoguanosine·isocytidine. (C) The base pair adenosine·pseudouridine that has a stabilizing backbone interaction. (D) Base pairs that result in significant helical distortions to the typical A-form RNA duplex: guanosine·uridine and inosine·uridine wobble pairs. Letter codes are shown above respective base pairs with hyphens representing Watson-Crick base pairing and middle dots representing wobble or modified base pairing. For simplicity, each base is shown capped with a methyl group where bases would connect to the sugar. The base pairs here are drawn to maximize base-base hydrogen bonding and to mimic the conformation of the canonical pairs. In solution within a duplex, it is likely that these pairs adopt this conformation; however, it is possible that they adopt a different hydrogen bonding pattern and/or slightly different conformation in varying sequence contexts.

Figure 2.

Process of generating nearest-neighbor geometries from MD trajectories. (A) An average structure is calculated from a 1 ns MD trajectory on each RNA duplex using only heavy (non-hydrogen) atoms. (B) Hydrogen atoms are added to the RNA duplex according to residue templates. Terminal base pairs and backbone are removed. (C) Bases are capped with H atoms at N9/N1 for purines/pyrimidines, respectively. From each MD simulation, two nearest neighbor geometries are generated for QM calculations.

Figure 3.

Schematic of bases involved for calculated hydrogen bonding energies and intra- and inter-strand stacking energies. E_HB was calculated for base interactions B1–B4 and B2–B3. Stacking energies (E_stack) include the two intrastrand (B1–B2 and B3–B4) and two interstrand base stacking combinations (B1–B3 and B2–B4).

Figure 4.

Experimental NN free energy parameters for Watson–Crick nearest neighbor combinations versus computational NN binding energies (ωB97X-D3 with CPCM (water, α = 1.5 Å)) from NN geometries obtained from MD simulations. Average fiber diffraction data was used to benchmark QM methods as described in Supplementary Figure S2. However, because our method to obtain modified base pair NN geometries comes from MD simulations, it was necessary to run the same MD protocol on all A–U and G–C base pair NN combinations to ensure consistency in energy derivation. Therefore, eight A–U and eight G–C computational NN free energies (Supplementary Tables S5 and S6) were mapped to experimental NN free energy parameters (19) using a simple linear regression.

Figure 5.

Predicted versus experimental nearest neighbor free energies. The dotted line has a slope of one and an intercept of zero, representing ideal correlation between predicted and experimental data. Watson–Crick (WC), inosine·cytidine (I·C), isoguanosine·isocytidine (iG·iC), 2,6-diaminopurineriboside·uridine (DAP·U) and inosine·uridine (I·U) experimental NN parameters were taken from references (19), (24), (25), (26) and (23), respectively. Individual panels for each modified base pair can be found in Supplementary Figure S4.