A Multidimensional B-Spline Correction for Accurate Modeling Sugar Puckering in QM/MM Simulations

Abstract

The computational efficiency of approximate quantum mechanical methods allows their use for the construction of multidimensional reaction free energy profiles. It has recently been demonstrated that quantum models based on the neglect of diatomic differential overlap (NNDO) approximation have difficulty modeling deoxyribose and ribose sugar ring puckers and thus limit their predictive value in the study of RNA and DNA systems. A method has been introduced in our previous work to improve the description of the sugar puckering conformational landscape that uses a multidimensional B-spline correction map (BMAP correction) for systems involving intrinsically coupled torsion angles. This method greatly improved the adiabatic potential energy surface profiles of DNA and RNA sugar rings relative to high-level ab initio methods even for highly problematic NDDO-based models. In the present work, a BMAP correction is developed, implemented, and tested in molecular dynamics simulations using the AM1/d-PhoT semiempirical Hamiltonian for biological phosphoryl transfer reactions. Results are presented for gas-phase adiabatic potential energy surfaces of RNA transesterification model reactions and condensed-phase QM/MM free energy surfaces for nonenzymatic and RNase A-catalyzed transesterification reactions. The results show that the BMAP correction is stable, efficient, and leads to improvement in both the potential energy and free energy profiles for the reactions studied, as compared with ab initio and experimental reference data. Exploration of the effect of the size of the quantum mechanical region indicates the best agreement with experimental reaction barriers occurs when the full CpA dinucleotide substrate is treated quantum mechanically with the sugar pucker correction.

Figure 1

Transition states of the RNA transesterification reaction models and RNase A system. (Only nonhydrogen “heavy” atoms are shown). (a) Minimal EthO···SP model used to construct 2D adiabatic PESs in the gas-phase; optimized MP2 transition state shown. (b) Full dinucleotide (CpA) nonenzymatic model. (c) RNase A enzyme complexed with CpA. Labeling of the O2′, O3′, and O5′ positions follows the conventions used in RNA, and lengths are given for the nucleophile (P–O2′) and leaving group (P–O5′) bonds.

Figure 2

QM regions used in the QM/MM simulations. TQM and FQM denote truncated and full QM regions, respectively. The histidines are RNase A residues and they are, therefore, not present in the nonenzymatic simulations of CpA.

Figure 3

2D gas-phase PESs of the EthO···SP nonenzymatic model transesterification reaction. The structures are the rate-limiting transition states. The coordinates are the nucleophile bond P–O2′ and the leaving group bond P–O5′. Stationary points along reaction paths are shown as black dots. Relative potential energies are calculated with respect to reactant energy minimum from each model (kcal/mol). Values in the parentheses are pseudorotation phase angle P_θ (deg) and puckering amplitude A_r (deg), respectively.

Figure 4

Pseudorotation cycle of the five-membered ring along with a furanose structure. The pseudorotation can be characterized by using the (A_r,P_θ) polar coordinates or (Z_x,Z_y) Cartesian coordinates. The relation between coordinate systems and proper torsions ν can be expressed as Z_x = (ν₁ + ν₃)/(2 cos(4π/5)), Z_y = (ν₁ − ν₃)/(2 sin(4π/5)), P_θ = arctan(Z_y/Z_x), and $A_r=\sqrt{Z_x^2+Z_y^2}$, where ν₁ and ν₃ are chosen to use the same starting phase angle as the convention ( $P_{\theta }{(}_2^3\mathrm{T})=0°$). The pseudorotation names are shown along the wheel for Exo (_nE), Endo (ⁿE), and Twist ( ${}_n{}^{m}T$, _nT^m, and ⁿT_m) conformations, where m and n denotes the O4′, C1′, C2′, C3′, and C4′ atoms in respective order; the superscript/subscript position stands for the atom above/below the ideal flat five-membered ring; the atom n on the left side of T has larger displacement from the ideal flat ring than the atom m on the right side.

Figure 5

1D PMF profiles of the nonenzymatic CpA model reaction. R5 and R2 denote the breaking bond P–O5′ and the forming bond P–O2′, respectively. Experimental values are indicated by dotted lines.

Figure 6

Coordinate overlay of approximate transition state structures from the nonenzymatic CpA model reaction simulations. The hydrogen atoms and cytosine and adenosine groups are removed to improve clarity.

Figure 7

Sampling, probability distributions, and statistical values of the P_θ and A_r (deg) pseudorotation parameters for the CpA cytidine sugar ring. MIN and TS denote reactant and transition states, respectively.

Figure 8

Sampling, probability distributions, and statistical analysis of the P_θ and A_r (deg) pseudorotation parameters from the RNase A reactant state C5 sugar ring.

Figure 9

2D PMFs of the enzymatic RNase A transphosphorylation reaction and their 1D projections. Stationary points along reaction paths are shown as black dots. R_GA and R_PT denote the reaction coordinates for general acid step (R_N119–H−R_O5′–H) and phosphate transfer step (R_P–O5′−R_P–O2′), respectively. LowR and UpL paths stand for the lower-right and upper-left asynchronous concerted mechanisms, respectively. RNaseA_TQM and RNaseA_FQM denote truncated and full CpA as QM regions in QM/MM simulations (Figure 2).

Figure 10

Sampling, mean values and standard deviations of the P_θ and A_r (deg) pseudorotation parameters from the C5 sugar ring in the RNase A reactant state (MIN) and transition states (TS_l and TS_u). TS_l and TS_u denote transition states along lower-right and upper-left asynchronous concerted paths, respectively (Figure 9). RNaseA_TQM and RNaseA_FQM mean truncated and full CpA as QM regions in QM/MM simulations (Figure 2).