|Home | About | Journals | Submit | Contact Us | Français|
Force field parameters specifically optimized for residues important in the study of RNA catalysis are derived from density-functional calculations in a fashion consistent with the CHARMM27 all-atom empirical force field. Parameters are presented for residues that model reactive RNA intermediates and transition state analogs, thio-substituted phosphates and phosphoranes, and bound Mg2+ and di-metal bridge complexes. Target data was generated via density-functional calculations at the B3LYP/6-311++G(3df,2p)//B3LYP/6-31++G(d,p) level. Partial atomic charges were initially derived from the CHelpG electrostatic potential fitting and subsequently adjusted to be consistent with the CHARMM27 charges and Lennard-Jones parameters were determined to reproduce interaction energies with water molecules. Bond, angle and torsion parameters were derived from the density-functional calculations and renormalized to maintain compatibility with the existing CHARMM27 parameters for standard residues. The extension of the CHARMM27 force field parameters for the non-standard biological residues presented here will have considerable use in simulations of ribozymes, including the study of freeze-trapped catalytic intermediates, metal ion binding and occupation, and thio effects.
Over the last several decades a wealth of data has accumulated that demonstrates the central role RNA catalysis plays in many biological processes. Starting in the late 1970’s, it was shown that RNA could catalyze fairly complex biological reactions in ribonuclease P1,2 and Tetrahymena3,4 with an efficiency that rivaled many protein enzymes. These discoveries sparked a wave of interest in the scientific community focused on unraveling the details of how RNA enzymes (ribozymes) function. An understanding of the catalytic mechanisms of ribozymes, and their relation to sequence and tertiary structure, is opening up a variety of new frontiers. In biomedical technology, gene expression inhibitors that target viral and genetic diseases5 such as HIV6 and cancer7 are being developed, and new biotechnologies such as RNA chips8 and allosteric molecular switches in nanodevices9 are being explored.
A common biological reaction catalyzed by several prototype ribozyme systems, such as the hammerhead,10,11 hairpin,12,13 and hepatitis delta virus14,15 ribozymes, involves cleavage of a phosphate group through a transesterification reaction16,17 (Figure 1). In this reaction, the 2′OH of the ribose sugar ring becomes activated via deprotonation and makes an in-line attack to the adjacent 3′-phosphate along the phosphodiester backbone. The attack produces a trigonal bipyramidal phosphorane intermediate/transition state that is accompanied by an inversion of configuration about the phosphorus center as the exocyclic P-O5′ bond is cleaved. The product of the transesterification is a 5′-OH terminus and a 2′,3′-cyclic phosphate. Additional information about this mechanism has been obtained via kinetic isotope studies18-20 and chemical modifications such as thio substitution21-23 at the scissile phosphate, although the mechanistic interpretation of these experimental studies remains a topic of discussion and some debate. Time-resolved x-ray crystallography10,24,25 has become a powerful tool to elucidate structural information at different stages along the catalytic reaction coordinate that provides valuable insight into ribozyme activity. However, this area is considerably challenging due to the difficulty of trapping a reactive intermediate and obtaining quality crystals, as well as uncertainties due to the nature of effects that arise from the crystallization conditions. For these reasons, there is considerable interest in the development of theoretical methods that can aid in the refinement and interpretation of existing experimental data, and provide structural insight into systems where such data is not yet available.
Molecular simulation, along with experimental structural data, provides an avenue for the characterization of ribozyme dynamics in solution and refinement of key mechanistic details. Molecular simulation force fields for nucleic acids continue to improve26-31 and a variety of simulations involving ribozymes have been carried out in recent years.32-40 In order to study the structure and dynamics of different catalytic states along the reaction path of a ribozyme, however, reliable empirical force fields parameters must be developed for the transition states and reactive intermediates of these reactions. Furthermore, in order to use molecular simulation to aid in the interpretation of experimentally measured thio effects, parameters for thio-substituted phosphate and phosphorane models and their interactions with metal ions are required. In the present work, new force field parameters for residues important to the study of RNA catalysis are derived from density-functional calculations to be consistent with the CHARMM2741 all-atom empirical force field. These parameters will allow molecular dynamics (MD) simulations of ribozymes in reactive states to be performed to study the structure and dynamics that lead to catalysis.
The first four summations are quadratic terms that give rise to energy penalties for geometrical deviations about equilibrium coordinate values. The variables b, θ, S, and are the bond length, bond angle, Urey-Bradley 1,3-distance, and improper torsion angle coordinates, respectively, while b0, θ0, S0, 0 and Kb, Kθ, KS, K are the corresponding force field parameters for the equilibrium geometries and force constants, respectively. The fifth summation is a trigonometric term that adjusts the periodic dihedral angle rotational barriers. The coordinate χ is the dihedral coordinate, n determines the periodicity, χ0 is a phase factor, and Kχ is the amplitude force constant. The terms involving a sum over atom pairs i, j < i (neglecting non-bonded exclusions) are the non-bonded van der Waals/Lennard-Jones (L-J) and electrostatic terms. The parameters i,j and R0,ij are the van der Waals well depth and minimum distance between the ij atom pair, respectively, and are by default calculated via the Lorentz-Berthelot combining rules45 from the corresponding 1-body parameters as and R0,ij = R0,i + R0,j, respectively, although this default can be explicitly overridden using the NBFIX (non-bond fix) option if fine tuning of specific pairwise interactions is desired. The last summation is the electrostatic energy determined between atomic partial charges qi, which normally are calculated with a unit dielectric constant, = 1, as in the present work.
Parameterization of the CHARMM27 force field was based on ab initio and experimental data for small molecules,42 as well as macromolecular simulation data of DNA and RNA46 in order to more accurately capture experimental condensed phase properties. This was an important improvement to the CHARMM22 force field47 that was parameterized more heavily on the basis of small molecule target data alone. While this previous approach succeeded in capturing the target data of the selected small compounds, later simulations of duplex DNA in solution showed some disagreement with experiment in regard to the relative conformational stability of A-DNA versus B-DNA.48,49
The CHARMM27 force field was developed using a self-consistent, step-wise strategy.31,42 Partial atomic charges for the CHARMM27 force field were initially obtained from a Mulliken population analysis of the HF/6-31G(d) wave function. These were then adjusted to better reproduce scaled HF/6-31G(d) TIP3P50 water interaction energies, experimental dipole moments, and heats of sublimation where available. Equilibrium geometry parameters for the selected small molecules were initially optimized to reproduce the experimental geometries of microwave, electron diffraction and/or x-ray crystal survey data where possible. Equilibrium geometries and force constants were optimized iteratively until satisfactory fitting to the target data was achieved. Once this optimization was complete, iterative adjustment of the charges and L-J parameters was coupled with the internal parameters until overall convergence was reached. A survey of RNA and DNA crystal structures from the Nucleic Acids Database (NDB)51 was taken so that fitting to target macromolecular experimental properties, such as sugar puckering phase and dihedral angle distributions, could be made through crystal MD simulation. During this stage, in the original CHARMM27 all-atom empirical force field parameterization for nucleic acids,42 dihedral angle parameters were adjusted to lower or raise energy barriers such that target condensed phase properties were better reproduced. While this sometimes sacrificed the quality of fitting to small molecule ab initio data, it accomplished the goal of more accurately reproducing experimental condensed phase properties. This macromolecular fitting was coupled with the internal parameter optimization until satisfactory convergence was achieved.
Due to the lack of high-resolution experimental data for RNA reactive intermediates and chemically modified nucleic acids, the parameterization of the present study is based solely on ab initio data, the optimized strucutures for which are shown in Figures Figures2,2, ,3,3, ,4,4, and and5.5. Parameter optimization for modified CHARMM27 RNA, Mg2+, and OH− residues was based on density functional theory (DFT) calculations for training sets of small molecules that represented the desired target systems. Density-functional calculations were performed using the Becke three-parameter hybrid functional combined with the Lee-Yang-Parr exchange-correlation functional (B3LYP)52,53 with the 6–31++G(d,p) basis set for geometry and frequency calculations followed by single-point electronic structure refinement with the 6–311++G(3df,2p) basis set in a manner analogous to recent studies of biological phosphates.54-58 The B3LYP model employed here neglects proper treatment of long-range dispersion interactions; however, these weak dispersion interactions are small in comparison with the highly polar (e.g., hydrogen bonded) and ionic interactions investigated here, for which B3LYP has been demonstrated previously to be generally reliable. All DFT calculations were performed using the GAUSSIAN0359 package. Partial atomic charges were based on the CHelpG60,61 method (CHarges from ELectrostatic Potentials using a Grid) which fits atomic centered point charges to the molecular electrostatic potential. The CHelpG charges were also constrained to reproduce molecular dipole moments. Similar approaches to force field charge determination have been applied and validated previously.60,62 These calculations were performed with the 6–311++G(3df,2p) basis set, and the CHelpG charges were subsequently modified to be consistent with the original force field by taking into account only the charge differences from standard CHARMM residues, and making further adjustments for polar and non-polar hydrogens and integer charged groups in accord with the original parameterization procedure31,42 (see Section 3.2 for complete details). The CHelpG charges used in the present work to obtain the charge differences with respect to similar standard CHARMM residues are different than the Mulliken charges that were used as a starting point for optimization in the original CHARMM27 force field, and in somees these charge modelsdiffer significantly. However, it should be emphasized that the Mulliken charges As in the original CHARMM27 force field, transferring charges from the small model compounds to the nucleic acid fragments was accomplished by adding the charge of the removed hydrogen atom to the heavy atom from which it was deleted. Once charges were obtained, Lennard-Jones (L-J) parameters were determined by fitting CHARMM residue-TIP3P water interaction energies and distances to the density-functional target values. This was accomplished through scaling of the interaction energies and shifting of the binding distances such that the DFT values matched those of the original CHARMM27 force field for unsubstituted dimethyl phosphate (DMP−). Although the original CHARMM27 force field parameterization fixed waters to the TIP3P geometry in ab initio calculations, it was found in this work that the relative differences in geometries predicted by the DFT level of theory used compared better with the HF/6-31G(d) used in the original force field when the waters remained unconstrained (see Table S2 in Supporting Information). Internal parameters were fit to the density-functional geometries and energies for the model compounds.
The CHelpG charges used in the present work to obtain the charge differences with respect to similar standard CHARMM residues are different than the Mulliken charges that were used as a starting point for optimization in the original CHARMM27 force field. In general, the CHelpG charges are larger in magnitude, and due to the constraints, preserve the molecular dipole moments. The Mulliken charges on the other hand are considerably basis set dependent, particularly when diffuse functions are used, and generally smaller in magnitude. It should be emphasized that the Mulliken charges used as a starting point for parameter optimization with CHARMM27 were often significantly altered upon optimization to obtain intermolecular interactions with TIP3P water molecules, the final charges in many cases being closer to the CHelpG charges (e.g., for the non-bridging phosphate oxygens).
RNA generally has more accessible “folded” conformations than DNA63 and often conformational deformation in the phosphate backbone is requisite for catalysis.64 The CHARMM27 force field contains dihedral parameters for unsubstituted RNA based on ab initio gas-phase calculations of a variety of test compounds and adjusted with respect to A and B form DNA and RNA within the Nucleic Acid Database (NDB).51 These adjustments included the lowering of torsional energy barriers in certain regions after the observation that direct parameterization to the ab initio torsion profiles resulted in simulated dihedral distributions that were inconsistent (e.g., too rigid) compared with those from a large NDB survey of RNA and DNA.42 In the case of nonbridging thio-substitutions, no parameters exist. Therefore, attention was given to the C-O-P-O (i.e., α and ζ in Figure 6) dihedral parameters because of their role along the phosphate backbone and the possibility of significant change upon thio-substitution. As in CHARMM27, dimethyl phosphate was used as the model compound for native RNA as well as singly and doubly thio-substituted dimethyl phosphate to determine the effects of the nonbridging sulfur(s) on this torsion. The dimethyl phosphate and thio substituted torsional potential energy surfaces were generated using B3LYP/6-311++G(3df,2p)//B3LYP/6-31++G(d,p) with all degrees of freedom relaxed except for the dihedral of interest and the symmetric ζ dihedral, which is fixed in the trans conformation to induce symmetry and facilitate parameterization (Fig. 7).
As in the original force field parameterization,42 a self-consistent step-wise optimization approach was taken that involved the iterative adjustments of L-J and internal parameters (not including torsion parameters) until convergence of the fitting function was obtained. Initial force constant values were taken from harmonic fitting to potential energy surface scans of bonds and angles with distortions of 0.2 Å and 2.0o respectively (See example in Figure 7). Starting geometries for the new CHARMM* model compounds were taken directly from the density-functional results and read into CHARMM to perform an Adopted Basis Newton-Raphson (ABNR) minimization65 until a gradient of <10−6 kcal/mol·Å was reached. The results of the CHARMM minimization with the current parameters were used to construct a χ2 fitting function. In the present work, this function took the form:
Where the , , and are the CHARMM energy-minimized bond lengths (Å), bond angles (degrees), torsion angles (degrees) and interaction energies (kcal/mol), respectively, and , , and are the corresponding DFT reference values. The σ parameters, have the same units as their associated geometrical or energetic quantities (indicated by superscript), the inverse squares for which serve as weights in the unitless χ2(λ) merit function. A non-linear optimization procedure is then used to minimize the χ2(λ) merit function with respect to the CHARMM* parameters (this work) denoted by the vector λ. The σ parameters (and hence the weights) were adjusted empirically over the course of a stepwise minimization procedure in order to drive the parameter search in the most realistic and robust region of parameter space (see Supporting Information for for more details). In most cases, the last step-wise minimization procedure used to arrive at the final set of parameters was a quadratically convergent direction set minimization method66 that does not require explicit gradient information. Sets of conjugate directions were generated with the algorithm due to Powell in a series of successive line minimizations achieved by parabolic interpolation. The convergence criterion was 10−4 on the relative change in the function value with respect to its minimum value after a series of line minimizations. While this method works very well for finding the local minima of a multi-variable function, it is limited in its ability to find the absolute minimum of such a function. For this reason, the minimization procedure was carried out at a large number of starting points for each variable until the best overall minima could be found. Once a minimum was reached the torsion angle parameters were adjusted empirically to obtain a reasonable balance between the trends in the relative ab initio curves (native versus sulfur substituted) while maintaining consistency with the standard CHARMM27 torsion profile which does not correspond directly to an ab initio curve. Course grained adjustment of the parameters (0.05 kcal) was seen to be sufficient to provide an acceptable balance. The torsion parameters were adjusted iteratively with the parameters of the non-linear optimization. The complete topology and parameter files of this work are available as Supporting Information.
RNA enzymes show diversity in the number as well as the role of metal ions they require. Metal ions can have both structural and chemical roles in catalysis. Possible (chemical) catalytic roles for metal ions include acting as:11,17
One way for a metal ion to act as a general base is to lose a proton from a coordinated water to form a metal hydroxide. This OH− has been postulated to play the role of abstracting a proton from a ribose 2′-OH,67 although the pKa difference between a phosphate coordinated Mg2+ and a ribose 2′OH remains a topic of current interest.68 In one postulated two metal ion mechanism for the hammerhead ribozyme (HR), it has been proposed that a metal hydroxide exists as a μ-bridging OH− between two magnesium ions.39 Other studies suggest a single metal ion mechanism with an outer-sphere coordination of the Mg2+ ion to the 2′OH69 and still others suggests a non-bridging two metal ion mechanism with inner-sphere binding of Mg2+.70 Clearly, MD and QM/MM simulation could help resolve these mechanisms and the catalytic role played by Mg2+ ions.
The force field parameters for Mg2+ and OH− are important to obtain accurate geometries and energies for Mg2+ ligand binding and the formation of di-metal Mg2+ complexes containing a bridging hydroxide. Current non-bonded parameters exist for Mg2+ in the CHARMM27 force field46 and OH− parameters are available from Lee et al.71 However, these parameters are not optimal for the specific problems that involve di-metal bridges or RNA phosphate binding. Hence, the parameters for OH− and Mg2+ have been re-optimized specifically for these types of systems.
DFT calculations have been run on several model systems55 which formed the training set for the Mg2+ and OH− parameterization. Each structure of the training set (Figure 3) contained one or two Mg2+ ions with a first coordination sphere combination of H2O, OH− and/or DMP− ligands. In the structures with two Mg2+ ions present, a OH− ion was placed in a μ-bridging position between the two metal ions. CHelpG charges were calculated for OH− and the L-J parameters for Mg2+ and OH− were iteratively adjusted until the best fit with geometry and water binding energy was obtained (Table 1). Mg2+ coordination distances were significantly improved with the new force field parameters for ligand binding with H2O, OH− and DMP−. This improvement was most pronounced for the di-metal bridging hydroxide structures which lengthen the Mg2+…OH− interaction by almost 0.30 Å and contract the Mg2+…OH…Mg2+ angle by almost 10 degrees. The new parameters also improved the Mg2+ water binding energy by approximately 1.0 kcal/mol. In the current study, only complexes that involve inner-sphere coordination of Mg2+ were used in the parameterization, and transferability to the important case of outer-sphere coordination remains to be tested with simulation.
The first step of transesterification for a RNA phosphate is deprotonation of the sugar 2′OH. This creates a nucleophile of sufficient strength to attack the phosphorous of an adjacent or nearby phosphate group.16 A “patch”42 has been prepared for the standard CHARMM27 RNA nucleotide which removes the 2′OH hydrogen and modifies the 2′O− charge, L-J interactions and geometry of nearby atoms. Parameterization of this patch is discussed below.
To derive new charges for a deprotonated nucleotide residue, a ribose methyl phosphate in its protonated and deprotonated states was optimized using DFT. It was found that an “in-line” geometry minimum did not exist in the gas phase without the solvated ribozyme environment. In order to parameterize a deprotonated ribose methyl phosphate (ribose MeP−) residue to a structure that models the in-line phosphate attack conformation, a minimization was performed with the corresponding α, and ζ dihedrals of the phosphate-sugar backbone fixed to the minimum energy geometry found for the phosphorane structure described in Section 3.3 (Figure 2). In this way, the backbone of the ribose MeP− was constrained to be “in-line”. CHelpG charges, constrained to reproduce molecular dipole moment for neutral molecules, were then calculated for the atoms of the optimized structures (Table 2). Since the original CHARMM27 charges (qo) were obtained via a slightly different procedure than those of the present DFT/CHelpG method (q), adjustments to the CHelpG charges (q’) were made to create charges more consistent with the constrained CHARMM27 force field parameterization. For instance, in the CHARMM27 force field, all non-polar hydrogens are constrained to have a charge of 0.09 e while polar hydroxyl hydrogens are assigned a charge of 0.43 e. The calculated CHelpG charges were adjusted for this by changing the value of the hydrogens to be consistent with the CHARMM27 force field while adding the difference created from this adjustment to the adjoining heavy atoms (Table 2). The differences between the adjusted CHelpG charges for the protonated and deprotonated ribose methyl phosphate were then calculated and used to modify the original CHARMM27 charges of a neutral nucleotide as follows:
where q’P and q’D are the adjusted DFT/CHelpG charges for the protonated and deprotonated structures, respectively, qo are the CHARMM27 charges, and are the new charges for the deprotonated structure. Final partial atomic charge parameters are shown in Table 2. It should be noted that constraints on the charge fitting to create unit charge groups within residues was not used in this method.
Once charges were calculated for the deprotonated ribose 2′O− oxygen, L-J parameters were determined by fitting to the 2′O− water coordination distance (1.713 Å) and interaction energy (-18.4 kcal/mol) obtained from DFT optimization of a ribose MeP− (Table 3). Bond and angle parameters were fit to the differences between the gas-phase DFT geometry of protonated and ribose MeP− (Table 4). Force constants were predicted by fitting to relaxed potential energy surface scans. All dihedrals except for H3′-C3′-C2′-O2′ were set to zero in the original protonated CHARMM27 structure and remained so in the new parameterization. The H3′-C3′-C2′-O2′, which has a Kχ of 0.195 kcal/mol/radian2 and 0.0° phase for the 3-fold term, was retained. Final fitting was performed for the residue patch within the RNA sequence UCA taken from the “early intermediate” active site geometry in the x-ray crystal structure of a hammerhead ribozyme.10 ABNR minimization was performed for this sequence with the U and A residues fixed and, where necessary to improve fitting, small adjustments to the internal parameters were made manually. The final results (Table 4) were able to reproduce the DFT structure and water binding energy very closely.
The deprotonated ribose 2′O− attack at the phosphate produces a phosphorane that is a transition state in the gas phase (Figure 2). Partial charges for the phosphorane atoms were determined similarly to those for the 2′O− ribose MeP− described above. However, due to the instability of dianionic phosphorane in the gas phase, DFT optimization was carried out with the axial P-O2′ bond length fixed to the P-O2′ bond length of a 2′O− ribose methyl phosphorane (ribose MePA2−) optimized in the aqueous phase (1.986 Å) using the Polarizable Continuum Model (PCM) solvation model.72,73 Comparison of similar neutral and anionic phosphorane structures optimized in both the gas and aqueous phases, suggests that the error for this constraint should be less than 0.02 Å for all but the axial P-O bonds which may have errors as high as 0.07 Å. Calculated CHelpG charges for hydrogen atoms were adjusted to 0.09 e with heavy atoms taking up the difference as before (Table 5). Final charges (q*TS) for the phosphorane structure were determined by adding the adjusted CHelpG charge difference (δq) between the ribose methyl phosphorane and the 2′O− ribose MeP− to the previously determined CHARMM* charges (q*D) for the the 2′O− ribose MeP− (see Table 2).
Bond length and angle equilibrium parameters for the ribose MePA2− were determined from fitting to the gas phase DFT optimization of the partially frozen dianionic ribose MePA2− structure described previously. Force constants for bonds and angles were calculated by fitting to relaxed potential energy surface (PES) scans for each P-O bond and O-P-O angle of a monoanionic phosphorane (which is stable in the gas phase). The final ribose MePA2− residue was fit within the RNA sequence UCA taken from the near in-line geometry in the x-ray crystal structure of a hammerhead ribozyme “late-intermediate”.74 ABNR minimization was performed for this sequence with the U and A residues fixed and where necessary, small adjustments to the internal parameters were made to improve fitting. Because the penta-coordinated geometry of the phosphorane is quite rigid, the dihedrals had negligible effect, and hence were set to zero. The final fitting results matched with the DFT results almost exactly (Tables 6).
Transesterification terminates in the creation of a 2′,3′-cyclic phosphate after exocyclic P-O5′ bond cleavage. The 2′,3′-cyclic phosphate is similar to a standard RNA ribose methyl phosphate residue except for the cyclization of the O2′-P-O3′-C3′-C2′ atoms (see Figure 2 and the parameter file in the Supporting Information). In the cyclic residue, the O2′ is bonded to the phosphorus in the same manner as the O3′. Due to the near symmetric equivalence between the O2′ and O3′ cyclic ribose phosphate sequences, the bond length, angle, and dihedral parameters are assumed to be identical. Therefore, O2′ was assigned the same atom type as O3′ (type ON2). Adjusted CHelpG charges for the 2′,3′-cyclic phosphate ring were calculated and force constants were determined through fitting to relaxed PES scans. Parameters were fit to the DFT data using ABNR minimization of the CHARMM* residue. New dihedrals were set to zero analogous those of the CHARMM27 ribose ring. Final parameters fit with very small errors with respect to the DFT results (Table 6).
To determine the role of divalent metal ions in a reaction, it is often useful to perform thio effect experiments.75,76 Changing the oxygen atom to a sulfur effects how strongly a divalent metal can bind to this position due to the differences between their sizes, polarizabilities, and bond lengths with phosphorous. While a “hard” ion like Mg2+ will bind very tightly to a correspondingly “hard” oxygen, it will bind much more weakly to a “softer”, more diffuse sulfur.77 If the divalent metal ion binding at a certain position is required for catalysis, a large decrease in the reaction rate should be seen upon thio substitution at this position. However, the interpretation of thio affect results are sometimes inconclusive21,23,78,79 especially when possible conformational changes induced by the sulfur are considered. Molecular simulation of thio-substituted phosphates and phosphoranes and their role on structure would aid in the interpretation of experimental thio effects.
A training set of phosphate compounds, bound and unbound to hexa-coordinated Mg2+, with both single and double sulfur substitutions was constructed (Figure 4). Adjusted CHelpG partial atomic charge calculations were calculated for both mono- and di-thio DMP− substitutions. Iterative optimization of the L-J parameters for sulfur was made until the best possible fit to the DFT calculated sulfur water-binding distance and relative water-binding energy between DMP− and di-thio substituted DMP− was obtained. Force constants and equilibrium geometry parameters were then optimized iteratively along with the L-J parameters until minimization of the χ2 function was achieved.
In order to improve the differential binding of Mg2+ to the non-bridging phosphoryl O and S atoms, a NBFIX term42 was used to parameterize specifically for Mg2+-S interactions; i.e., a specific value for the 2-body non-bonded van der Waals parameters was included explicitly for Mg2+-S rather than using the Lorentz-Berthelot combining rules to derive the 2-body parameters from 1-body Mg2+ and S parameters. The resulting fit was further improved by the introduction of three new atom types. Two for the non-bridging oxygen (ONS) and sulfur (SO) of a mono-substituted phosphate and one for the non-bridging sulfur (SS) of a di-thiosubstituted phosphate. The DFT training set average binding distance for a Mg2+-S was 2.55 Å as compared to the CHARMM* average value of 2.06 Å (Table 6). The shortened distances arose out of the need to strengthen the binding energy results that, in the absence of explicit polarization on the soft sulfur atoms, is considerably underestimated with non-polarizable force fields. In the case of water binding to the substituted sulfurs, the distances are in better agreement and the relative binding energies are all within 1.0 kcal/mol of the DFT values.
The ζ ab initio torsional potential energy surfaces for the native and nonbridging thio-substituted DMP− are shown in Figure 7 (top frame). The minima identified agree well with fully relaxed structures in previous work by Florian et al.80 A significant difference between the thio-substituted and unsubstituted surfaces is the minima at the staggered position, 180°, which is negligible for the native DMP−, but of almost equal energy to the gauche minima in the case of dithio-substituted DMP−. Also noteworthy is the shift of the gauche minima (70°) and gauche-staggered barrier with thio substitution. New torsional parameters, summarized in Table 6, are defined along the S-P-O-C and O-P-O-C dihedrals to qualitatively reproduce the changes indicated by the ab initio calculations, while maintaining consistency with the CHARMM27 force field.
A scan of the DMP− and thio-substituted torsions with the same constraints as the ab initio calculations using the CHARMM27 parameters and those developed in this work are shown in Figure 7 (bottom frame). As indicated in Foloppe et al.,42 the barrier between gauche-gauche states in the unsubstituted DMP− has been significantly lowered to reproduce results from experiment, while the gauche-staggered barrier is left relatively unchanged. The thio-substituted structures retain the lowered gauche-gauche barrier and the gauche minima is shifted as indicated by the ab initio calculations. The staggered conformation has been lowered in energy for both thio-substituted molecules. For the singly substituted case, the minima at 300° was also lowered in addition to the minima at 60°, but the relative energy of the two gauche and eclipsed conformations were maintained. The lowering staggered conformation may be of particular importance to catalytic differences in native and thio-substituted structures since conformational deformation is a prerequisite to reaction.32,63,81,82
L-J parameters for the non-bridging atoms of the penta-coordinated di-anionic phosphorane structures were parameterized using a training set which contained methyl ribose phosphorane (ribose MePA2+), mono-thio substituted ribose MePA-so2+ and di-thio substituted ribose MePA-ss2+, with and without water bound (see Figure 5). In the case of the unsubstituted phosphorane, it was required to freeze the P-O2′ bond as discussed in Section 3.3. Stable gas-phase sulfur substituted phosphoranes were used for the other structures. Both the non-bridging S and O atoms were reparameterized in this case due to the large differences in geometry and charge distribution of a phosphate and phosphorane. As before, equilibrium geometry values were fit to DFT structures with force constants for all modified bonds and angles determined through fitting to relaxed PES scans. Geometry fitting results are shown in Table 6.
Molecular simulation force field parameters are presented for a series of non-standard residues important in RNA catalysis. Parameters are based on density-functional calculations and developed to be consistent with the CHARMM27 all-atom empirical force field for nucleic acids. Parameters have been developed for an activated (2′O deprotonated) ribose phosphate representing an early reactive state, a penta-coordinate phosphorane intermediate/transition state model, and a 2′,3′-cyclic phosphate transesterification product. In addition, parameters for thio-substituted analogs important in the study of experimental thio effects, and specific single and di-metal Mg2+ complexes and μ-bridging OH− ion parameters have been developed. These parameters, which are optimized specifically for the present systems, will allow the simulation of reactive intermediates and experimentally modified residues important in the study of ribozyme mechanisms. Further validation and testing of these parameters in simulations is important, although this will likely be a lengthy and tedious endeavor due to the lack of available experimental structural data for the reactive intermediates and thio-substituted nucleic acids. Ultimately, it is the hope that simulations of these systems, together with experiment, will help paint a more detailed picture of the catalytic mechanisms of RNA catalysis.
DY is grateful for financial support provided by the National Institutes of Health (Grant GM62248), and the Army High Performance Computing Research Center (AHPCRC) under the auspices of the Department of the Army, Army Research Laboratory (ARL) under Cooperative Agreement number DAAD19-01-2-0014. Computational resources were provided by the Minnesota Supercomputing Institute. ADM is grateful for financial support provided by the National Institutes of Health (Grant GM51501).