Xabier Lopez

Annick Dejaegere

Fabrice Leclerc

Darrin M York

Martin Karplus

March

Nucleophilic Attack on Phosphate Diesters: A Density Functional Study of In-Line Reactivity in Dianionic, Monoanionic, and Neutral Systems

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Phosphate diesters play a fundamental role in biology, including their role as the backbone of DNA and RNA [1,2]. The chemical properties and reactivity of phosphates determine how these biomolecules are formed and cleaved, and therefore, phosphate diesters have been the subject of numerous theoretical and experimental studies [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26].

A first approximation to understand phosphate diester reactivity is to characterize the energetics of gas-phase model reactions for which a full quantum mechanical treatment is possible. Of particular interest are studies that cover nucleophilic attack of water or methanol on ethylene phosphate (EP), a model for the transphosphorylation and hydrolysis of RNA chains.

Previous computational studies on non-enzymatic phosphate hydrolysis reactions [3,9,5,11,27,12,13,17,23,24] have focused mainly on the so-called dianionic and monoanionic reaction mechanisms, in which either a negatively charged nucleophile (XO -; X=H,CH 3 ) or a neutral one (XOH; X=H,CH 3 ) attacks an unprotonated phosphate diester molecule (i.e., dimethyl phosphate or ethylene phosphate). Very few of these studies has addressed the effect of the sugar ring on the reaction. The pK a values of phosphates (typically below 3) suggest that they are ionized in aqueous solution around neutral pH, and hence are the most likely reactant species in nonenzymatic hydrolysis. However, the protonation state of the phosphate esters in enzymatic hydrolysis, especially of the phosphorane transition states and intermediates, are not clear. Recent experimental and theoretical results suggest that the phosphoranes exhibit significantly elevated pK 1 a (e.g., a value of 7.9 has been suggested for ethylene phosphorane [27] and 8.6 for P(OH) 5 [28]). Consequently, it is important also to characterize neutral reaction mechanisms since the charge state of the phosphorus species can play an important role in the enzyme catalyzed reactions.

In this paper, we extend the scope of previous computational work [16,29] to consider neutral reaction mechanisms and to characterize the effect of the sugar ring at a high level of theory and basis set. We have studied by means of density-functional theory the following set of reactions:

XO(H) + P O 2 (H)(-O -Y -O) → XO -P O 2 (H) -(O -Y -O(H)) (1)
where X can be either hydrogen or methyl, and Y is either C 2 H 4 or a ribose ring. These reactions are illustrated schematically in Fig. 1. We restrict the focus to in-line attack mechanisms, and do not consider mechanisms that involve pseudorotation. Although exploration and characterization of other mechanisms is important, the formation of the pentacovalent intermediate that results from in-line attack is generally the rate limiting step [30,31]. Calculations were performed using density-functional theory (B3LYP functional) at high basis set levels (6-31++G(d,p)/6-311++G(3df,2p)). Solvation was approximated by the linear isotropic PCM dielectric continuum model [32,33,34].. Non-enzymatic phosphoryl transfer reactions with quantum chemical methods have been studied extensively with quantum chemical methods in order to quantify the influ-ence of factors that regulate reactivity. Numerous studies have been performed in the gas-phase [4,3,5,17,35], and in solution [9,10,6,7,8,16,27,18,23,24,26,36,15,37]. These studies utilize different levels of electronic structure theory, models for solvation and other computational protocols that reduce the degree to which meaningful cross-comparisons between studies can be made. The present paper presents calculations for a series of reactions carried out at a consistent level of theory and treatment of solvation that helps unite and extend the results presented previously by different authors with different methods. [4, 5, 6, 7, 8, 9, 10, 11, 16, 27, 38, 39, 26, 23, 24, 40, ?] Consequently, the results should be of general utility and in particular provide a benchmark dataset for the optimization of semiempirical quantum models that can be used in the application of hybrid quantum mechanical/molecular mechanical methods to study enzymatic phosphate hydrolysis reactions in biology [29].

Methods

Gas-phase calculations

All the structures in the gas-phase were optimized with Kohn-Sham density-functional theory (DFT) methods, using the hybrid exchange functional of Becke [41] and correlation functional of Lee, Yang and Parr [42]. Frequency calculations were performed to verify the nature of the stationary points on the potential energy surface (PES);

i.e., that there were no imaginary frequencies for minima, and only one imaginary frequency for transition states. Geometries and frequencies were determined with a 6-31G++(d,p) basis set. Following geometry optimization, the electronic energies were refined using a larger 6-311++G(3df,2p) basis set. This basis set is similar to that used in the G2 method [43,44], which typically yield atomization energies, ionization energies and relative energies within 1 kcal/mol of the experimental values. All the calculations were done with the GAUSSIAN98 suite of programs [45].

The B3LYP/6-31G++(d,p) frequencies were used to evaluate the zero-point vibrational energy (ZPVE), thermal vibrational contributions to the enthalpy, entropy and

Gibbs free energy at T=298 K within the harmonic oscillator approximation in the canonical ensemble [START_REF] Hehre | Ab inito Molecular Orbital Theory[END_REF][START_REF] Cramer | Essentials of Computational Chemistry: Theories and Models[END_REF]. The rotational and translational energies were treated classically as (1/2)RT per degree of freedom. The thermodynamics relations and energy contributions calculated in the gas-phase are:

G = H -T S (2)
H = E + RT (3)
E = E 0 + E vib + E rot + E trans (4)
E 0 = E e + ZP V E (5)
where G, E, H, S and T are the Gibbs free energy, the energy, the enthalpy, the entropy and the temperature (T=298 K), respectively, and E e , E vib , E rot and E trans are the electronic, thermal vibrational, rotational and translational energy components, respectively. All quantities except E e and ZP V E have an explicit temperature dependence.

The most serious deficiencies of the evaluation of gas-phase free energies in this work has to be attributed to dianions, for which sufficiently flexible diffuse basis functions are required for accurate calculations. On the other hand, the magnitude of the gas-phase entropy calculated as the sum of the translational, rotational and vibrational entropies calculated using the ideal-gas, rigid rotor and harmonic oscillator approximations often leads to a large overestimation of the entropy term for reactions in solution.

This is so since the harmonic approximation underestimates the entropy contribution from low frequency modes. Besides, the consideration of infinitely separated reactants with fully free translational and rotational degrees of freedom as the reference state overestimates the entropy penalty associated to the approach of two reactants in solution, since there is an overestimation of loss of translational and rotational degrees of freedom upon approach of the two reactants. Finally, one should also consider that there is a neglect of the changes in solute entropy upon the transfer from the gas-phase to solution.

Solvation calculations

Solvent effects were treated using the polarizable continuum model (PCM) [32,33,34]. In this model, the solute molecule is embedded in a cavity of unit dielectric that is surrounded by solvent modeled by a dielectric continuum. The cavity was defined using the UAHF sets of radii [START_REF] Barone | [END_REF]. The solute and solvent are coupled by a reaction potential of the dielectric medium in response to the solute charge distribution.

The polarization of the solvent is represented by a charge density σ introduced on the surface S of the cavity surrounding the solute, and the corresponding reaction field potential takes the form

φ(r) = S d 2 s σ(s) |r -s| (6)
To check the sensitivity of solvation corrections to the method and basis set choice, the solvation calculations were performed with two different electronic structure methods (HF and B3LYP) and two different basis sets (6-31G(d) and 6-31+G(d)), for a total of four combinations of theory/basis set (HF/6-31G(d), HF/6-31+G(d), B3LYP/6-31G(d) and B3LYP/6-31+G(d)). The free energy in solution was calculated as a solvation free energy correction to the gas-phase free energy as:

G aq = G gas + ∆G sol (7)
where ∆G sol is the solvation free energy calculated as the the difference between the electronic energies in solution and the gas-phase, plus the self-energy of the solvent reaction-field along with an empirical term that account for the non-electrostatic contributions to the solvation free energy (sometimes referred to as the "cavitation" term) [32,49]. The solution free energies of the gas-phase stationary points from Equation 7are then used to estimate the relative energies of the various stationary points in solution, ∆G i aq with respect to the energies in solution of the separate reactants. The relative values in solution with respect to the reactant are less affected by the basis set and method choice than absolute solvation energy values, so that the shape of the potential energy surface in solution is similar for the different methods tested. Therefore, and for the sake of brevity we chose only one method , HF/6-31+G*, to report aqueous free energies; the results obtained with the other methods are given as supplementary material.

Due to the large number of structures considered, the solvation free energies were calculated at the gas-phase B3LYP/6-31G++(d,p) geometries rather than performing geometry optimizations in solution.

It is known that geometries can be affected by inclusion of solvent effects [16] and a complete treatment of solvation should include geometry optimizations and consideration of key explicit solvent molecules that form structural hydrogen bonds.

Nonetheless, the gas-phase geometry/continuum dielectric approximation applied here can provide important qualitative insight into how solvent affects the energetics of the different reactions [9,10,11]. Another drawback of the present calculations is the lack of explicit solvent water molecules as effective catalyst for the intramolecular proton transfers from phosphoryl oxygens to departing nucleophiles/leaving groups in monoanionic and neutral reactions. Notice that water-mediated proton transfer has been claimed to facilitate formation/cleavage of phosphoranes [50].

Nomenclature

To label the different structures along the reaction profile, we introduce the following notation. There are reactants (REAC), hydrogen bonded complexes (HB), products (PROD), and pentacovalent phosphorane intermediates (INT) that are identified as minima on the potential energy surface, as well as several transition states (TS1, TS2 and TSINT). The name TS1 refers to the transition state for the approach of the nucleophile to the phosphate reactant, while TS2 refers to the transition state for the endocyclic cleavage which connects the phosphorane intermediates to the products.

TSINT refers to transition states that connect different phosphorane intermediates (see Result section for details).

The number of stationary point identified on the potential energy surface differs for the dianionic, monoanionic and neutral reactions. In the case of dianionic reactions, there are only REAC, TS1 and PROD structures. In the case of monoanionic and neutral reactions there are more possibilities and the nomenclature is detailed below.

Monoanionic Reactions

In the case of water attack on ethylene phosphate, we have considered two type of intermediates (Int a and Int b ), which differ in the orientation of the equatorial proton (See Fig. 4). Considering that the five membered ring in ethylene phosphate is not planar, the two equatorial phosphoryl oxygen are not equivalent, so that there are two Int a and two Int b structures. However, preliminary calculations showed almost equivalent energetics (within 1 kcal/mol) and similar structures when protonating one or the other oxygen, so that only one structure will be discussed here.

For the reaction with ribose ethylene phosphate, there are more isomers, according to whether the protonated phosphoryl oxygen is exo or endo with respect to the ribose ring. These two possibilities were specifically considered, and we distinguish both structures with an exo and endo superscript (See Fig. 4).

Neutral Reactions

In the neutral reactions, the stationary points identified are also reactants (REAC), hydrogen bonded complexes (HB), transition states (TS), phosphorane intermediates (INT) and products (PROD). However, the number of different isomers is larger, as the equatorial oxygens are both protonated and each proton can assume different orientations. In this case we have introduced two subscripts which take into account the orientation of each of the equatorial protons. The symbol d refers to "down" orientation of the equatorial proton according to the structures of Fig. 5, meaning that this proton is oriented towards the oxygen of the nucleophile. The subscript u refers to the "up" orientation, i.e., proton oriented towards the axial endocyclic oxygen (i.e. the oxygen of the leaving group). As shown in Fig. 5, we use two subscripts in each case, the first one refers to the orientation of the proton bound to the exo phosphoryl oxygen, and the second one corresponds to the equatorial proton bound to the endo oxygen (exo or endo with respect to the ribose ring).

Results and Discussion

The results of DFT calculations for the series of phosphate hydrolysis reactions are presented and discussed. The section is organized into three main subsections that focus on dianionic, monoanionic and neutral reactions, respectively. Each of these subsections examines structures, and thermodynamic quantities in the gas-phase and in solution. Particular attention is paid to the geometries of the phosphorane interme-diates and the transition states corresponding to the formation of the exocyclic P-O apical bond, the overall rate limiting transition state in all cases, and to the specific protonation states that give rise to possible phosphorane isomers.

Dianionic reaction mechanism

Gas phase profile

The dianionic reaction mechanism corresponds to the attack of the anionic reactant (hydroxide or methoxide) on the negatively charged phosphodiesters, ethylene phosphate and 2',3'-ribose ethylene phosphate . These are the small molecule models for mechanisms frequently proposed for enzymatic phosphate hydrolysis, where the attack of an activated nucleophile (for example, deprotonated hydroxyl) occurs on the phosphate diester in its physiological anionic protonation state. This reaction has been studied previously on simpler model systems lacking the ribose ring [3,9,10,11,16,18] and the global features of the reaction profiles presented here (cf Figure 3 and Table 1) are in accord with previous studies in that hydrolysis proceeds via direct in-line nucleophilic displacement with inversion of configuration at the phosphorus. The rate limiting transition state occurs early in the reaction, with distances between the phosphorus and the incoming nucleophile oxygen around 2.645 Å (cf. Table 2) and an endocyclic P-O ax distance of 1.846 Å. The reaction proceeds directly from the transition state to the products without formation of stable pentacovalent dianionic phosphorane intermediates. In the P-O ax distance range where stable intermediates are observed for the monoanionic and neutral reactions (see below), the potential energy surface is rather flat, but no stationary point was identified. Thus even if shallow intermediates exist, they are not kinetically significant [51]. This type of mechanism is consistent with kinetic 18 O isotope effects [50] for the cleavage reaction of Uridine 3 ′ -m-nitrobenzyl phosphate (i.e., the reverse transesterification reaction to the methanolysis of cyclic phosphate diesters discussed in the present work), supportive of a concerted type mechanism but with a late transition state (that is, equivalent to our early-type transition state for the methanolysis reaction of the cyclic phosphodiester). In addition, on the basis of leaving group effects, highly negative values for β lg (slope of Bronsted-type plots)

in the case of the hydroxide ion catalyzed transesterification cleavage of uridine 3'alkylphosphates [52,53,54] have been reported, suggesting that at the transition state the departure of the corresponding alkoxide group is almost complete, which is again in agreement with an early transition state for the methanolysis of cyclic phosphodiesters. Quite interestingly, a break in the Bronsted-type free energy relationship has been recently reported [54] for the alkaline cyclisation of uridine 3'-phosphate esters at pK a values of the leaving group (the nucleophile in our case) lower than the pK a of the 2' oxyanion nucleophile. This convex break in the Bronsted plot has been taken as an evidence of a transient intermediate.

The dianionic reaction has a very high activation energy barrier in the gas-phase, it is of the order of 81 (ethylene phosphate) to 87 (2',3'-ribose ethylene phosphate

) kcal/mol (see Table 1). The free energy barrier is higher by 9 to 10 kcal/mol, which reflects the loss of translational and rotational degrees of freedom in the di-anionic complexes, the inclusion of vibrational entropic contributions does not significantly affect the shape of the reaction path. The high activation barrier results largely from the coulombic repulsion between the two negatively charged ions in the gas-phase [9,10,11]. The presence of the ribose ring lowers the activation barrier by 5 kcal/mol with respect to ethylene phosphate, presumably due to its better ability to delocalize the electronic charge and reduce the coulomb repulsion of the dianion.

Overall, the hydrolysis reaction is endothermic in the gas-phase. The nature of the attacking group (methoxide vs. hydroxide) influences the reaction energetics for the hydrolysis; i.e., when methoxide is the nucleophile, the reaction is more endothermic by over 10 kcal/mol (see Table 1). This can be linked to the greater stability of the methoxide reactant relative to hydroxide. The effect of the ribose ring is fairly small with respect to the barrier height and overall reaction energetics. It can be noted that the backward reaction, i.e., that from the products to reactants 2',3'-ribose ethylene phosphate and methoxide is a model for the transphosphorylation step of RNA hydrolysis. From the data presented here, the activation energy for this particular reaction is 32.8 kcal/mol, which is significantly less than that for the hydrolysis. The activation energy for the backward hydrolysis reaction is higher (42.4 kcal/mol) when hydroxide is the leaving group, as expected, since hydroxide is a worst leaving group than methoxide in bearing a negative charge in the gas-phase.

Solution profile

To obtain a first approximation to the reaction profile in solution, we estimated the solvation free energies (G sol ) of the gas-phase stationary points with the continuum dielectric (PCM-UAHF) methods as discussed in section 2.2 (Tables 1,S9, S10 and S11).

It can be seen that for the dianionic reaction, the solvation free energy of the separate reactants is significantly less favorable than that of the transition state and product.

This is in accord with the expectation that a dianionic species is more favorably solvated than two monoanionic reactants [11]. The large favorable solvation of dianionic transition states leads to a significant reduction of the activation barrier in solution. The barrier is lowered to 45-55 kcal/mol depending on the nature of the reactants. In the gas-phase, the activation barrier for the 2',3'-ribose ethylene phosphate reaction was somewhat lower, while as the transition state for the 2',3'-ribose ethylene phosphate compounds are bulkier, they are less well solvated and therefore the estimated barriers in solution for the 2',3'-ribose ethylene phosphate reactions are similar, or even somewhat higher than those obtained for ethylene phosphate. A detailed comparison of the solution barriers for the various reactants and nucleophiles would necessitate obtaining stationary points in solution. Nevertheless, the current data indicate that the barriers are similar for all reactions.

These data also highlights the importance of the consideration of the solvent to understand the hydrolysis reactions, as noted in earlier studies. [10,11]. For the reverse (i.e., transphosphorylation) reaction, the effect of solvation on the activation barrier is quantitatively much less important than for the hydrolysis step, as the reactants are also dianionic; solvation reduces the activation barrier by about 10 kcal/mol. The difference in solvation free energies with the different methods employed tend to cancel out when calculating relative energies. This can be seen in Fig. 7, where we depict the values of the solvation free energies for the different species calculated at HF/6-31+G*, B3LYP/6-31G* and B3LYP/6-31+G* versus the HF/6-31G* values. The inclusion of diffuse functions tend to slightly increase the solvation free energy (in absolute value), whereas the B3LYP method has a larger effect and gives smaller (in absolute value) solvation free energies than Hartree-Fock in all cases. However, these trends are consistent for all the structures treated, and hence, the barriers calculated in solution with the different methods are similar. As stated in the method section, only HF/6-31+G* numbers are included in Tables 1, 3 and 6.

Monoanionic reaction mechanism

Gas phase profile

Monoanionic reactions correspond to either the attack of an activated nucleophile (i.e.

OH-or CH3O-) on the neutral phosphate or to the attack of a neutral nucleophile (water or methanol) on the phosphate in its usual charge state (-1). Monoanionic reactions have been proposed in enzymatic systems [START_REF] Breslow | Proc. Natl. Acad. Sci[END_REF] where the attack of an incoming charged nucleophile would be accompanied by a proton transfer from a nearby general acid to the phosphate, which would facilitate the attack.

For the reactions considered here in the gas-phase, identical stationary point are obtained for the attack of an anionic nucleophile on the neutral phosphate and the attack of the neutral nucleophile on the anionic phosphate, so that except in the reactant stage, the potential energy surfaces are identical. The initial step in the gas-phase reaction for the attack of water (or methanol) on the monoanionic phosphates is the formation of an H-bonded complex between the reactants (see Table 3 and Fig. 4).

The stability of the hydrogen bonded complex is quite significant in the gas-phase due to favorable electrostatic interactions (-15.5 kcal/mol in energy) but the entropic contribution to complex formation is unfavorable, as expected, and hence in free energy, the stability is only -3.3 kcal/mol (Table 3). The geometry of the complex differs if water or methanol is used as nucleophile, with water two slightly elongated hydrogen bonds are formed with the non-bridging phosphoryl oxygens (approximately 2.1 Å), whereas CH 3 OH forms a single strong hydrogen bond (approximately 1.75 Å). If the reactants used for initial gas-phase geometry optimization are the neutral phosphate and anionic nucleophile, no stable complex is formed where the phosphate oxygen is the H-bond donor. This is expected from the relative gas-phase proton affinities of the reactants; i.e. when the neutral phosphate and anionic nucleophile interact, a proton transfer occurs and the only stable H-bonded complex that can be identified correspond to the neutral nucleophile donating the H-bond to the anionic phosphate oxygen.

As the reaction proceeds to the rate-limiting transition state, the proton forming the hydrogen bond is transferred from the nucleophile to one of the non-bridging phospho- In this sense, it has been shown that in the case of phosphorane monoanionic triesters, for which intramolecular proton transfer from the phosphoryl oxygen to the leaving group is not an option, isomerization is 10 5 times faster than the cleavage.

However, in phosphorane diesters isomerization is only 1 to 2 orders of magnitude faster than cleavage. Besides, also consistent with a proton transfer to the leaving group concerted with the rate-limiting bond rupture, the experimental [START_REF] Kosonen | [END_REF] β lg of the cleavage is only moderately negative, -0.59, as compared with the high negative value for the hydroxide-ion-catalyzed reaction of phosphodiester monoanions [52], -1.28, in which the alkoxide departs unprotonated.

In general, there are large departures from the ideal 180 o in the O ax -P -O ax angle in the transition state; values between 159 o and 164 o are obtained (Table 4).

Regarding the distance between the nucleophile and the phosphorus in TS1, they are all very similar, 2.235 -2.265 Å with slightly longer P-OX distances for methanol than for water. For the attack on 2',3' ribose ethylene phosphate, the P-OX distances in the transition state are shorter than for ethylene phosphate. The transition states TS1's of the monoanionic reactions are more associative than the transition states of the dianionic reactions, in which the values of P-OX and P-OC distances at TS1 were longer than in the present case. Due to the geometry of the phosphate reactant the two non-bridging equatorial oxygens are not equivalent, so that two transition states are identified (labeled exo and endo for 2',3'-ribose ethylene phosphate as a reactant, see Figure 4). However, there is very little difference in the activation barriers and overall geometries between the two structures.

The gas-phase energy barrier between the reactants (taken as the neutral nucleophile and charged phosphate) and rate limiting transition state are 21-25 kcal/mol, while the corresponding free energy barriers are 32-35 kcal/mol (see Table 3). The barriers with respect to the stable hydrogen bonded complex are of course higher, but in free energy the difference is only about 3 kcal/mol. For instance, in the case of a methanol attack on 2',3' ribose ethylene phosphate the gas-phase free energy barrier with respect to separated methanol and charged phosphate is 34.2 kcal/mol for T S1 Exo , while the barrier with respect to the hydrogen bonded complex is 37.8 kcal/mol. It is also possible to consider taking the anionic nucleophile and neutral phosphate as the reference energy. However, these reactants are much higher in energy than the charged phosphate and neutral nucleophile (by 50-65 kcal/mol in the gas-phase, see Table 3). In solution, the difference in energy is less (See Table 3, the anionic nucleophile and neutral phosphate are only about 10-15 kcal/mol higher than the negative phosphate and neutral nucleophile) ) but the lowest energy species remains the neutral nucleophile and anionic phosphate, as expected. It is thus clear that in the gas-phase, the activation barrier is significantly lower for the monoanionic reaction than for the dianionic one.

The potential energy surface for the monoanionic reaction is more complex than for the dianionic one. Indeed, as the reaction proceeds from the exocyclic transition state, stable pentacovalent compounds are formed (i.e., phosphoranes). They are the intermediates of the reaction and they have a distorted trigonal bipyramid structure.

The central phosphorus atom is bound to 5 oxygens, two of them occupy axial positions (O ax ) and other three occupy equatorial positions. These intermediates are higher in energy than the reactants by 13-15 kcal/mol (see Table 3). For the monoanionic reaction of 2',3'-ribose ethylene phosphate four of these intermediates can be identified (cf. Figure 4) that differ in the position and relative orientation of the proton bound to the equatorial oxygens. Indeed, the proton can be bound to each equatorial oxygen, as these are not equivalent (labeled exo and endo, cf. Figure 4). The proton can then be oriented either towards the incoming nucleophile oxygen (Int a ) or towards the departing leaving group oxygen (Int b , Figure 4).

These stable intermediates occur at different stages on the reaction pathway. In the transition state of nucleophilic attack (TS1), the proton that has been transferred from the incoming nucleophile is still oriented towards it, and this orientation is kept in the subsequent intermediate (Int a ). A rotation of the equatorial proton then occurs leading to the transition state (TSINT, cf. Figure 4) where the proton is situated in the equatorial plane. This transition state is associated with a significant barrier of about 10 kcal/mol with respect to the intermediates, and show similar gas-phase energies to those of the transition states for nucleophilic attack (see Table 3). This rotation of the proton bound to the equatorial oxygen is an interesting feature of the gas-phase potential energy surface. The barrier for rotation of the proton is wide enough that the probability of tunneling is negligible.

In the case of water as the nucleophile, the rotation of the proton is coupled to pseudorotation of the phosphorane intermediate. This type of coupling between the rotation of the P-OH bonds and pseudorotation has also been identified in the past [58].

Therefore, a transition state for pure proton rotation cannot be identified as a stationary point on the potential energy surface. A discussion of pseudorotation mechanisms [25] is beyond the scope of the present paper, and we did not characterize this transition state further. A similar feature had been noticed by Tole and Lim [5,6] in their early studies of hydroxyl ion attack on methyl ethylene phosphate, a different but related monoanionic system. There was found to occur from the transition state that has an axial hydroxyl group to place it in an equatorial position.

In the case of methanol as a nucleophile a transition state that corresponds clearly to the rotation of the proton around the P-Oeq bond is identified. The attack of methanol on ethylene phosphate was previously studied by Lim and Tole, at the HF/3-21+G* level and a transition state for proton rotation was observed by them [59].

The proton rotation leads to a different intermediate (Int b , Figure 4) where the proton is now oriented towards the leaving group oxygen, thus facilitating its departure through the final (TS2) transition state. The intermediates subsequent to proton rotation Int b are identified both in the case of OH-and CH3O-as nucleophiles.

The final (TS2) transition states correspond to the breaking of the axial endocyclic P-O bond (see Fig. 4). The bond breakage is accompanied by a proton transfer to the leaving group oxygen. The proton is closer to the equatorial oxygen and further from the leaving group axial oxygen in TS2 than in TS1 for exocyclic cleavage. (see Tables 4 and5). In fact the distance between the O phos. and the hydrogen is only slightly elongated with respect to what is observed in the stable intermediate. Correspondingly, the distance between the hydrogen and the leaving apical oxygen is significantly longer As for exocyclic attack, two TS2 are identified (exo and endo, cf. Figure 4), they are similar in energy and they are lower in energy than the other TS by 10 kcal/mol, and thus not rate limiting (cf. Table 3 and Figure 4).

From this transition state, the rupture of the endocyclic P-O axial bond occurs, together with the transfer of the hydrogen on the leaving group oxygen. This yields a product with a charged phosphate group and neutral ribose ring, as expected. In the gas-phase, an intramolecular H bond remains formed in this product between the ribose hydroxyl and a phosphoryl oxygen. This internal hydrogen bond stabilizes the reaction product in the gas-phase, so that overall, and contrary to the dianionic case, the reaction in exothermic. More extended geometries in which the internal hydrogen bond is broken are higher in energy by around 10 kcal/mol (data not shown).

Solution profile

As for the dianionic case, some important changes are observed when ∆G i aq and ∆G i gas are compared, which highlight the consideration of the solvent to understand these reactions. For the monoanionic species, the free energy of solvation of the gasphase stationary points on the potential energy surface is generally less favorable than the solvation of the separate reactants (see Table 3), as expected since the complexes have the same charge as the reactants and have larger radii.

The consequence is that the activation barrier for the monoanionic reactions increases by about 5 to 10 kcal/mol when passing from gas-phase to solution. The magnitude of this solvation effect is smaller than was observed for the dianionic species (see Tables 3 andS12). Therefore, the main effect of solvation is to bring the acti-vation free energies for the monoanionic and dianionic reactions much closer than in the gas-phase. Thus, for example for CH 3 O(H) attack on 2',3' ribose ethylene phosphate, we have a gas-phase free energy barrier of 92.98 kcal/mol and 34.18 kcal/mol for dianionic and monoanionic reactions, respectively. However, in solution, the barriers are 45.05 kcal/mol and 43.25 kcal/mol with solvent contributions evaluated at B3LYP/6-31+G* level of theory. Thus the barriers differ by about 60 kcal/mol in the gas-phase but only about 2 kcal/mol in solution. These results are in agreement with previous ab initio studies [3,9,5,11,27,12,31].

Solvation data presented here are single point energetic estimates on gas-phase structures, and explicit geometry optimization in solution could result in somewhat different structures. Region of the gas-phase potential energy surface where structures are likely to be affected by solvation are the stable H bonded complex forming at the beginning of the gas-phase reaction, which is unlikely to form in solution, as there will be competition from water molecules forming H-bonds with both reactants. The transition state for proton rotation that connects the two intermediates is better solvated than the other transition states which is expected in view of its smaller size and hydrogen atom pointing towards the solvent and so the rotation profile could be affected by solvation. (cf. Figure 4 and Tables 3 andS12).

Neutral reaction mechanism

Gas phase profile

The neutral reaction corresponds to the attack of water (or methanol) on the protonated phosphate. Although phosphates will be anionic at physiological (neutral) pH, acid catalyzed hydrolysis can proceed through the neutral mechanism at very low pH (i.e., for instance the pK a of Dimethyl Phosphate is 1.29 [27]). Moreover, phosphorane intermediates are more basic than phosphates [27,19], so that neutral phosphoranes can be relevant to reactions at near neutral pH. Relative energies with respect to the separate neutral reactants (protonated phosphate and water/methanol) can be found in Table 6, relevant geometrical information are found in Tables 7 and8. Full information on solvation free energies and aqueous free energy barriers with various theoretical levels can be found in Tables S14 andS15 in the supplementary material.

The reaction profile for the neutral reaction resembles in many aspects that of the monoanionic reaction mechanism. However, as both the nucleophile and phosphate are protonated, the number of isomers for the transition states are twice that observed for the monoanionic reactions, since the two protons that become bound to the equatorial oxygens can assume two different orientations (see Figure 5 and discussion below). The first step of the neutral reaction in the gas-phase is the formation of a stable hydrogen-bonded (HB) complex. However, this kind of complexes are kinetically irrelevant, specially in aqueous solution where positive free aqueous energies are calculated with respect to the infinitely separated reactants. The reaction then proceeds to the TS for exocyclic cleavage. Four distinct TS are identified when a ribose ring is present. All of them correspond to a nucleophilic attack on the phosphorus atom coupled with a proton transfer from the nucleophile to one of the equatorial oxygens. This type of mechanism in which P-O bond formation/rupture is coupled with proton transfer between phosphoryl oxygen and nucleophile has also been considered based on the β lg values of leaving groups on the cleavage of RNA neutral phosphodiesters [53],

which are only negative and thus, suggest that the leaving group in these reaction departs as an alcohol. In the case of a ribose ring present, we classified the two equatorial oxygens as being exo or endo with respect to the ribose ring, and therefore, there are two exo and two endo transition states (see Fig. 5), classified according to the equatorial oxygen that receives the proton in the transition state. For each of the exo and endo transition states, there are two possibilities, depending on the orientation of the other equatorial proton (the one that is already present in the phosphate reactant). As a result, in two of the TS1 transition states, T S1 exo dd and T S1 endo dd , both protons attached to the equatorial oxygens points towards the nucleophile, while in the two other transition states, T S1 exo du and T S1 endo ud , one proton points towards the nucleophile and the other towards the leaving group. The energies of the different transition states with respect to the reactants (see Table 6) do not differ much, even when comparing the reaction with different nucleophiles and with/without ribose ring. All the gas-phase free energy barrier lie within 36-38 kcal/mol. Larger differences are observed in the geometries (see Table 7 and8). When both equatorial protons are pointing towards the nucleophile, the P-Oaxial distance to the nucleophile O is longer, around 2.1 Å , than when one proton points up and the other down, around 1.95 Å (cf. Table 7). These P-O axial distances are shorter than in the monoanionic (longer than 2.2 Å see Table . 4)

or dianionic case (around 2.6 Å , Table 2). An important conclusion, therefore, is that the degree of associative/dissociative character is highly influenced by the protonation state of the system, augmenting the associative character of the transition states as the total charge decreases. The geometry for proton transfer in the TS also differs between the monoanionic and neutral cases. In the neutral TS1's, the proton from the incoming nucleophile is almost equidistant from the nucleophile oxygen and the phosphoryl oxygen (see Table 8), while in the monoanionic case (Table 5), the distance to the equatorial O is much shorter than to the nucleophile O. It is remarkable to note that β lg values for the transesterification cleavage of neutral phosphodiesters [53] are less negative, namely -0.12, than for monoanionic phosphodiesters, -0.59, consistent with a higher degree of proton sharing between the nucleophile and the phosphoryl oxygen at the neutral transition states than in monoanionic reactions. The reaction then proceeds to the formation of stable neutral phosphorane intermediates. There are four distinct neutral phosphorane intermediates depending on the orientation of the protons at the phosphoryl equatorial oxygens (see Figure 5), labeled as INT du , INT dd , INT du and INT uu . We notice however, that barriers for the connection of these intermediates are small and much lower in every case than the barriers implied by TS1 or TS2 (exo and endocyclic cleavage) type transition states. All the phosphorane intermediates lie very close in energy. In general the following trend is observed: i) alternant orientation of the protons, ud or du, is favored ii) intermediates resultant from water attack show somewhat lower relative energies than intermediates resultant from methanol attack, specially if for the intermediates that show the equatorial protons oriented towards the nucleophile ("down" orientation). In general our data agree with the results obtained in [27]. All these intermediates can be formally connected by transition states corresponding to different single equatorial proton rotations as shown schematically in Fig. 6. However, we found that in many cases these rotations were coupled with pseudorotation of the whole structure at the transition state (TSINT transition states). Transition states with a significant pseudorotation character show O ax -P -O ax angles around 160 o . A full discussion of these types of transition states are beyond the scope of the present paper, an interesting discussion of pseudorotation processes in oxyphosphoranes can be found in ref [25]. Barriers for the connection of these intermediates are small and much lower in every case than the barriers implied by in-line TS1 or TS2

(exo and endocyclic cleavage) type transition states. From the neutral intermediates, the reaction can proceed to endocyclic cleavage, only if at least one proton is oriented towards the leaving group so as to protonate it. Four transition states that have the proton in this orientation are identified (cf. Figure 6). Again, when a ribose ring is presented, the transition states can be classified exo or endo, depending on whether the proton transfer to the leaving group is exo or endo. Overall, these TS for endocyclic cleavage are lower in energy than those for exocyclic cleavage (by about 6 kcal/mol) and thus they are not rate limiting. In the neutral transition state for endocyclic cleavage, the proton is shared between the leaving group oxygen and one equatorial oxygen, as was already observed for the exocyclic TS (cf. Table 8 and Fig. 5). The products are finally formed by cleavage of the P-OC bond and simultaneous transfer of one of the phosphoryl hydrogen to the leaving group. An internal H-bond is maintained, as for the monoanionic case. The reactions energetics of the neutral and monoanionic reactions are similar in the gas-phase (see Tables 3 and6).

Solution profile

Solvation has a relatively minor effect on the reaction energetics of the neutral reaction (see Table 6), as compared to monoanionic or dianionic reactions. Solvation, nevertheless, tend to favor the separated reactants, and therefore, we obtain higher relative energies in solution than in the gas-phase. The hydrogen bonded complexes formed at the beginning and end of the gas-phase reaction are significantly destabilized by solvation, suggesting that they would not occur in solution; a geometry optimization would be needed confirm they remain present in water. As for the transition states and intermediates, they also show higher relative energies in solution, but the favorable solvation free energy of separated reactants is somewhat balanced by the higher polarization of the P-O bonds in phosphorane structures than in tetravalent phosphate species.

Overall the shape of the neutral reactions are very similar to that of the monoanionic reactions. However, the energies and free energies of the neutral phosphorane intermediates with respect to the neutral reactants differs from the monoanionic case (see Table 6), the neutral phosphoranes being significantly more stable, between 14-18 kcal/mol for neutrals, and between 25-28 kcal/mol for monoanionic ones. Therefore, an important consequence is that the barrier for the dissociation of the intermediates into products/reactants are significantly higher for the neutral than for the monoanionic reaction path. This along with the low barriers for the interconversion among intermediates, suggest that neutral phosphorane intermediates will show sufficient lifetime to pseudorotate and therefore the endocyclic cleavage reaction can be proceed from various possible intermediates.

Reaction pathways of cyclic phosphodiesters

The hydrolysis of phosphodiesters is generally thought to follow an associative mechanism (see [60,61,62] and ref. cited therein) and the stationary points identified here for the hydrolysis reactions of cyclic phosphodiesters clearly agree with that view. Significant differences in the structures of the transition states according to the total charge are nevertheless observed, and it is of interest to examine the degree of bond breaking/forming at the different stationary points along the reaction pathway. We compared the bond lengths between the nucleophile/leaving group and phosphorus in the transition states and intermediates with their values in the product or reactants, respectively (see Fig. 8). Significant bond formation occurs between the nucleophile and phosphorus in the rate limiting transition states for hydrolysis, while the distance to the leaving group is only slightly increased with respect to the distance in the reactants. This is the hallmark of an associative mechanism [63], confirming early experimental studies in which it was remarked that dissociative mechanism is highly unlikely for the hydrolysis of phosphodiesters [60]. More recent 18 O kinetic isotope studies [50] and determination of β lg Bronsted-type slopes for the cleavage of hydroxide-ion-catalyzed RNA phosphodiester bonds [52,53,54,64] also support an associative mechanism even for the dianionic reaction. However, the degree of bond forming to the nucleophile varies significantly with the total charge, with the dianionic reaction having the longest P-OX distance, and hence the least associative transition state. For the monoanionic and neutral reactions, the transition state occurs with more bond formation to the nucleophile than in the dianionic reaction. It can, however, be seen that the elongation of the bond lengths in the transition states with respect to the reactant/product value (cf.

Fig. 9a) remains small in all cases (it is at most 30% for the dianionic case), which is in agreement with the associative character of the reactions. Interestingly, for the neutral reactions, the orientation of the equatorial protons (cf. Figures 4 and5) influences the geometry of the transition states: when both protons are oriented towards the nucleophile, the transition state occurs earlier (i.e. with less bond formation between the nucleophile and phosphorus) than when one of the protons is oriented towards the leaving group. This is in accord with the Hammond postulate [65,16]. For the neutral and monoanionic reactions, the reaction proceeds through the formation of phosphorane intermediates, and in those cases, the bonds to both the nucleophile and leaving group are only at most 10% longer than in the reactants/products. The formation of these stable phosphoranes is also a clear indication of an associative mechanism. For the neutral and monoanionic reactions, departure of the leaving group occurs from the intermediates, while the bond to the nucleophile remains close to the product value.

The elongation in the transition states for leaving group departures (TS2, cf. Figures 4 and5) is at most 25% of the reactant/product state. The link between the geometries and energies is presented in Figure 9. The calculated activation barriers with respect to the separate reactions are similar for all reactions. The intermediate clearly appear as low energy structures and increased stability of the neutral phosphoranes is seen. For the neutral and monoanionic reactions, the rate determining transition state is always that for attack of the nucleophile rather than departure of the leaving group. For the dianionic mechanism, the transition state also corresponds to attack of the nucleophile rather than departure of the leaving group. This seems to be coherent with experimental kinetic analysis [54,64] of isomerization and cleavage reactions of RNA phosphodiesters. Quantitative estimates of reaction barriers would require more detailed calculations, including geometry optimization in solution and possibility of water-mediated proton transfer, which has been shown to play an important role facilitating the proton transfer in monoanionic and neutral reactions [50].

Concluding Discussion

We have present an extensive computational study of a variety of cyclic phosphodiester hydrolysis reactions. The calculations and discussion focused on the most relevant stationary points of each potential energy surface to explain the in-line mechanism for this important set of biological reactions. From our systematic studies, we confirm that the transition state corresponding to the initial attack of the nucleophile on the phosphate is the rate-limiting transition state of the reaction. This trend is preserved regardless of the particular nucleophile, the presence/absence of a ribose ring in the model phosphate or the inclusion of solvent effects. The total charge of the system is the main factor that determines the energetic trends in the gas-phase: the largest barriers are for the dianionic mechanism, while the monoanionic and neutral reactions have similar energetics and smaller barriers. Another important effect of the total charge is the existence and stability of phosphorane intermediates. No gas-phase intermediates are observed for the dianionic reaction. Monoanionic intermediates are found in the gas-phase, but the barrier for their dissociation to the products is very low. Neutral phosphorane intermediates, on the other hand, have significant barriers for dissociation to both reactants and products, which suggest that they could have a significant lifetime in the gas-phase. Phosphoranes that have kinetically significant lifetimes may undergo protonation/deprotonation from solvent, or pseudorotation that may alter the stereochemical outcome of the reaction and lead to products involving migration [2,31,1]. The data presented here are thus coherent with the observation that isomerization reactions leading to migration products are catalyzed under acidic conditions [1] (which favor neutral phosphoranes). The total charge also influences the geometry of the transition state and intermediates.

All the reactions presented here correspond to associative mechanism, yet differences are observed with the total charge. The rate limiting transition state becomes gradually more associative, with P-OX bonds decreasing by 0.69 Åand P-O ax by 0.17 Å when going from the -2 to the neutral reaction. Concurrently, the O ax -P-O ax angles approach linearity as the negative charge of the system decreases. Thus, the average O ax -P-O ax for monoanionic intermediates is 161.9 o , whereas for neutral intermediates is 178.4 o . Solvent has a dramatic effect on the relative energies with respect to the infinitely separated reactants and serves to equalize the barriers for the dianionic, monoanionic and neutral mechanisms. Our calculations also highlight similarities with the hydrolysis mechanism of non-cyclic phosphate such as monomethyl [36,15] and dimethyl phosphate [66]. Indeed, in an early study of the mechanism for associative monomethyl phosphate hydrolysis in gas-phase and in solution [36,15], it was observed that the monoanionic and neutral reaction mechanism proceed by proton transfer from the nucleophile to the phosphate reactant, and that solvation raises the activation barrier for the neutral and monoanionic reactions, while it leads to a significant lowering of activation free energy for the dianionic case. A recent study[66] of the dianionic and monoanionic mechanism for the hydrolysis of dimethyl phosphate reached similar conclusions. Our work also showed that for the monoanionic and neutral reactions, proton transfer from the nucleophile to the phosphate is an essential step of the reaction. Moreover, we characterized various isomers of the monoanionic and neutral transition states and showed that for ribose ethylene phosphate, both exo/endo proton transfer lead to very similar barriers.

We also showed that isomerization of the protonated phosphoranes play an essential role in the reaction pathway and are associated with significant energy barriers in the gas-phase. This is a point that has important mechanistic implications for enzyme catalyzed reactions i.e. it shows that for neutral and monoanionic reactions, the orientation of protons attached to the equatorial oxygen has an important role in facilitating attack/departure of the nucleophiles, and that the change in proton orientation is an component of the reaction path that should not be overlooked.

This work has implications for mechanistic interpretation of phosphate hydrolysis in biological system, particularly in RNA hydrolysis where cyclic phosphate intermediates are commonly encountered. In enzymatic systems, it is often postulated that a general acid protonates the phosphate and/or a base activates the nucleophile. Here it is seen that in the simple gas-phase model system constituted of phosphate and nucleophile, the nucleophile acts as the general acid and the phosphate as the general base.

These data provide important insight into the intrinsic reactivity of the phosphates in the absence of surrounding enzyme or solvent and therefore provide a reference point for understanding more complex systems in which catalysis occurs. It can also be taken as a high-level density functional theory dataset of reactions for which more approximate models, such as focused semi-empirical treatments, can be developed and tested. Table 5: Proton transfer between oxygens at the transition states. H -O X is the distance in ( Å) between the proton and the oxygen of the nucleophile (water or methanol), H -O P refers to the distance between the proton and the corresponding equatorial phosphoryl oxygen, H -O -CH 2 -is the distance between the transfered proton and the axial O -CH 2 -group. O-O is the distance between the oxygen proton donor and acceptor and O-H-O angle is the proton transfer in degrees, respectively.

H -O X H -O P H -O -CH 2 -O-O O -H -O H -O X H -O P H -O -CH 2 -O-O O -H -O HOH + EP- CH3OH + EP- T S1 1

  ryl oxygens. This proton transfer is almost complete in the transition state, as the distance between the incoming O and the proton (1.4 -1.5 Å) is significantly longer than the distance between the phosphoryl O phos. and the hydrogen (1.0 -1.1 Å). The H-O-H angle however deviates significantly from linearity (134-137 o , see 5). and the O --O distance (proton donor and acceptor) is between 2.308 and 2.325 Å . Experimental results on cleavage and isomerization of RNA phosphodiester bonds[START_REF] Kosonen | [END_REF]57,53], suggests a similar mechanism, in which the rate-limiting transition state corresponds to exocyclic cleavage of the phosphorane concerted with proton transfer to the departing alkoxide.

( 1 .

 1 70 to 1.75 Å) than in the exocyclic TS (1.4 to 1.5 Å). The H-O-H angles deviate even more from linearity in the endocyclic TS (120 o -125 o ), and the O-O distance between proton donor and acceptor are slightly longer, 0.1 Å . On the other hand, in the presence of a ribose ring, both exo and endo transition states are very similar, but the exo transition states are located slightly later in the reaction coordinate, being the P-OC distances at the transition state slightly longer for exo transition states, and correspondingly, the proton transfer is slightly more advanced at exo transition states than at endo ones.

Figure 1 :

 1 Figure 1: Scheme for the in-line mechanism in phosphate ester hydrolysis.

Figure 2 :

 2 Figure 2: Phosphate reactants: Ethylene phosphate (EP-), protonated ethylene phosphate (EPH), 2',3' ribose cyclic phosphate (REP), and protonated 2',3' ribose cyclic phosphate (REPH).

Figure 3 :

 3 Figure 3: Dianionic Potential Energy Surface. B3LYP/6-31++G(d,p) optimized geometries of the separate reactants, transition state and product for the in-line methoxide attack on ribose ethylene phosphate.

Figure 4 :

 4 Figure 4: Monoanionic Potential Energy Surface. B3LYP/6-31++G(d,p) optimized geometries of transition states and intermediates for the in-line methanol attack on ribose ethylene phosphate.

Figure 5 :

 5 Figure 5: Neutral Potential Energy Surface. B3LYP/6-31++G(d,p) optimized geometries of transition states and intermediates that corresponds to the in-line mechanism of the nucleophilic attack of methanol on 2p,3p-cyclic ribose ethylene phosphate. The dd, du , ud and uu superscripts refers to the orientation (up or down) of the equatorial protons.

Figure 6 :Figure 7 :

 67 Figure 6: Transition States connecting the intermediates in the neutral potential energy surface.

Figure 8 :

 8 Figure 8: Axial P-O interatomic distances with respect to reference values in the product (P) and reactant (R) at each TS1, INT and TS2 structures.

Table 1 :

 1 Dianionic reactions: Relative energies in the gas-phase with respect to the separate reactants for the XO -+ EP -/REP -reactions (X=CH 3 , H; REP=2',3'-cyclic ribose ethylene phosphate, EP -=ethylene phosphate). Numbers are in kcal/mol. E e corresponds to the electronic energy calculated at our best level of theory (see Methods). Contributions from Zero-Point Vibrational Energy, enthalpy and entropy are estimated at B3LYP/6-31++G(d,p). Last column for each reaction type corresponds to the relative free energy in aqueous solution ∆G aq , obtained from the sum of ∆G gas and the ∆∆G sol solvent contributions (See Methods) evaluated at HF/6-31+G* TS1 86.71 87.14 86.43 95.94 43.46 81.32 81.80 81.14 90.56 54.11 PROD 35.88 37.00 36.60 44.89 20.45 38.89 39.77 39.34 47.85 23.10 .31 87.13 98.33 44.07 81.97 82.12 82.05 92.98 45.83 PROD 44.10 46.34 46.51 56.64 19.83 49.16 51.14 51.33 61.63 25.18

		∆E e ∆E 0	∆H	∆G ∆G aq ∆E e ∆E 0	∆H	∆G ∆G aq
		Reaction: CH3O -+ EP -	Reaction: CH3O -+ REP -
	TS1	87.17 87

Reaction: HO -+ EP -

Reaction: HO -+ REP -

Table 2 :

 2 P-O axial bond distances ( Å) and angle (in degrees) for all the stationary points characterized for the XO -+ (R)EP -reactions using B3LYP/6-31+G(d,p) level of theory. P-OX P-OC O ax -P-O ax P-OX P-OC O ax -P-O ax

		HO-+ EP-	HO-+ REP-
	TS1	2.627 1.834	163.07	2.645 1.846	162.9
	PROD 1.700	91.4	1.700 4.264	103.9
		CH3O-+ EP-	CH3O-+ REP-
	TS1	2.627 1.834	163.07	2.645 1.846	162.9
	PROD 1.700 5.136	89.6	1.707 4.599	108.6

Table 3 :

 3 Monoanionic reaction: Relative energies in the gas-phase with respect to the separate reactants for the XOH + EP -/REP -reactions (X=CH 3 , H). Numbers are in kcal/mol. E e corresponds to the electronic energy calculated at our best level of theory (see Methods). Contributions from Zero-Point Vibrational Energy, enthalpy and entropy are estimated at B3LYP/6-31++G(d,p). Last column for each reaction type corresponds to the relative free energy in aqueous solution (∆G aq ), obtained from the sum of ∆G gas and the ∆∆G sol solvent contributions (See Methods) evaluated at HF/6-31+G*.

		∆E e	∆E 0	∆H	∆G ∆G aq ∆E e	∆E 0	∆H	∆G ∆G aq
			Reaction: HOH + EP -			Reaction: CH3OH + EP -
	HB	-15.27 -12.82 -13.36 -4.59 8.93 -14.65 -13.49 -13.22 -4.86 7.50
	T S1	22.88 23.73 22.23 34.26 42.14 21.84 21.29 20.63 32.96 41.94
	Int a	13.34 16.45 15.07 27.12 26.68 13.73 15.44 14.77 27.85 30.72
	T SINT	-	-	-	-	-	22.65 23.26 22.52 35.58 35.77
	Int b	13.76 16.29 15.24 26.32 25.88 13.73 15.44 14.77 27.85 30.72
	T S2	13.70 15.83 14.35 26.37 29.50 13.78 14.78 13.85 27.29 31.90
	PROD	-19.36 -16.33 -17.31 -6.76 0.88 -19.35 -17.52 -17.84 -6.32 6.72
	XO -+ EP H	61.02 60.98 60.96 62.06 13.71 52.11 49.71 49.72 50.24 15.61
			Reaction: HOH + REP -			Reaction: CH3OH + REP -
	HB	-14.51 -11.94 -12.53 -3.29 10.27 -14.00 -12.68 -12.47 -3.62 9.18
	T S1 Exo	24.48 25.07 23.66 35.33 39.53 23.64 22.80 22.24 34.18 43.96
	T S1 Endo	24.28 25.08 23.64 35.40 39.56 23.42 22.86 22.25 34.45 43.18
	Int Exo a	13.64 16.57 15.31 26.98 26.48 14.70 16.05 15.57 27.96 29.27
	Int Endo a	13.65 16.51 15.33 26.71 27.34 13.49 15.04 14.51 27.10 29.78
	T SINT Exo	-	-	-	-	-	23.24 23.60 23.01 35.69 34.78
	T SINT Endo	14.51 16.92 15.49 27.34 26.75 22.41 22.90 22.28 34.90 35.58
	Int Exo b	13.02 15.86 14.74 26.01 28.33 13.98 15.54 15.03 27.61 31.45
	Int Endo b	13.65 16.51 15.33 26.71 27.34 13.49 15.04 14.51 27.10 29.78
	T S2 Exo	13.99 15.93 14.56 26.19 30.51 14.23 14.96 14.23 26.86 30.96
	T S2 Endo	13.69 15.65 14.27 25.97 28.45 13.79 14.58 13.81 26.75 31.73
	PROD	-14.49 -11.93 -12.75 -2.77 8.27 -14.43 -12.96 -13.15 -2.04 9.55
	XO -+ REP H 65.79 65.45 65.00 66.98 14.15 56.88 54.18 53.76 55.17 16.06

Table 6 :

 6 Neutral reaction. Relative energies in the gas-phase with respect to the separate reactants for the XOH + EPH/REPH reactions (X=CH 3 , H). Numbers are in kcal/mol. E e corresponds to the electronic energy calculated at our best level of theory (see Methods). Contributions from Zero-Point Vibrational Energy, enthalpy and entropy are estimated at B3LYP/6-31++G(d,p). Last column for each reaction type corresponds to the relative free energy in aqueous solution (∆G aq ), obtained from the sum of ∆G gas and the ∆∆G sol solvent contributions (See Methods) evaluated at HF/6-31+G*. CH 3 OH + EPH HB -12.06 -9.40 -10.24 -0.55 10.81 -12.46 -10.60 -10.79 -0.62 12.20 T S1 dd 28.02 28.01 26.61 38.08 41.62 26.56 25.46 24.73 37.18 43.33 T S1 ud 25.64 26.45 24.72 36.92 42.60 25.33 24.96 23.99 37.12 43.71 S2 du 19.75 20.26 18.67 30.68 36.56 22.18 21.41 20.59 33.45 39.63 T S2 uu 15.10 15.64 14.03 26.00 37.95 13.71 13.55 12.57 25.66 39.90 PROD -9.48 -7.21 -7.82 1.35 9.80 -11.05 -9.42 -9.48 0.82 11.23 CH 3 OH + REPH HB -12.09 -9.46 -10.30 -0.66 8.87 -12.47 -10.62 -10.23 -2.43 8.25 T S1 Exo du 25.98 26.86 25.66 36.64 39.43 25.71 25.37 24.95 36.81 40.09 T S1 Exo dd 28.49 28.41 27.61 37.60 38.56 26.93 25.87 25.68 36.76 42.06 T S1 Endo dd 28.24 28.14 27.34 37.23 38.06 26.66 25.54 25.38 36.40 40.84 T S1 Endo ud 25.75 26.68 25.46 36.39 39.50 25.46 25.09 24.69 36.37 41.50 .00 19.99 30.56 34.44 19.53 19.15 18.79 30.41 36.38 T S2 Exo uu 15.68 16.20 15.20 25.69 34.62 14.01 13.89 13.48 25.24 36.93 T S2 Endo du 20.89 21.33 20.31 30.95 33.98 23.30 22.36 22.17 33.47 37.48 T S2 Endo uu 16.17 16.53 15.54 26.02 34.95 14.61 14.23 13.91 25.28 36.12 PROD -7.02 -4.27 -4.52 3.91 9.19 -8.68 -6.72 -6.32 3.30 11.22

			.430	1.072		2.325	136.7	1.467	1.050		2.322	134.5
		T S2		0.986	1.746	2.335	120.4		0.985	1.748	2.336	120.3
					HOH + REP-					CH3OH + REP-
		T S1 Exo	1.401	1.085		2.314	136.8	1.429	1.065		2.308	135.1
		T S1 Endo	1.428	1.073		2.325	136.2	1.465	1.052		2.323	134.1
	58	T S2 Exo T S2 Endo		0.992 0.986	1.699 1.730	2.408 2.408	124.9 122.5		0.992 0.985	1.695 1.737	2.410 2.411	125.5 122.2

Table 7 :

 7 P-O axial bond distances ( Å) and angle (in degrees) for all the stationary points characterized for the XOH + (R)EPH reactions using B3LYP/6-31+G(d,p) level of theory. P-OX P-OC O ax -P-O ax P-OX P-OC O ax -P-O ax

		HOH + EPH	CH 3 OH + EPH
	HB	3.247 1.626	135.2	3.240 1.627	127.6
	T S1 dd	2.137 1.637	173.2	2.091 1.641	172.6
	T S1 ud	1.956 1.678	179.2	1.936 1.685	178.4
	Int du	1.677 1.724	197.6	1.673 1.729	175.2
	Int dd	1.726 1.680	176.8	1.723 1.686	179.7
	Int ud	1.677 1.724	179.6	1.674 1.727	179.6
	Int uu	1.644 1.790	177.7	1.638 1.798	178.2
	T SInt u	1.641 1.734	157.8	1.656 1.754	172.7
	T SInt d	1.687 1.675	157.8	1.677 1.688	156.0
	T S2 du	1.641 2.070	175.6	1.626 2.130	172.5
	T S2 uu	1.600 2.289	173.9	1.589 2.320	173.9
	PROD	1.612 3.414	150.8	1.602 3.478	146.8
		HOH + REPH	CH3OH + REPH
	HB	3.253 1.624	132.2	.244 1.625	132.2
	T S1 Exo du	1.950 1.678	178.8	1.930 1.685	178.6
	T S1 Exo dd	2.141 1.636	173.5	2.099 1.642	173.2
	T S1 Endo dd	2.132 1.635	172.8	2.091 1.641	172.6
	T S1 Endo ud	1.952 1.679	179.4	1.931 1.686	178.5
	Int du	1.675 1.724	1.781	1.670 1.731	174.5
	Int dd	1.724 1.681	177.0	1.721 1.687	179.2
	Int ud	1.676 1.727	178.9	1.672 1.732	179.6
	Int uu	1.642 1.791	176.3	1.636 1.802	177.1
	T SInt Exo u	1.642 1.736	158.6	1.642 1.747	160.1
	T SInt Exo d	1.668 1.694	157.5	1.689 1.691	161.7
	T SInt Endo u	1.643 1.738	159.8	1.653 1.759	170.7
	T SInt Endo d	1.683 1.673	156.2	1.688 1.689	160.6
	T S2 Exo ud			1.629 2.093	178.2
	T S2 Exo uu	1.597 2.323	174.1	1.590 2.335	176.1
	T S2 Endo du	1.638 2.081	176.2	1.623 2.157	172.4
	T S2 Endo uu	1.597 2.323	174.1	1.586 2.359	174.1
	PROD	1.602 3.320	124.4	1.591 3.326	124.6

Table 8 :

 8 Proton transfer between oxygens at the transition states. H -O X is the distance in ( Å) between the proton and the oxygen of the nucleophile (water or methanol), H -O P refers to the distance between the proton and the corresponding equatorial phosphoryl oxygen, H -O -CH 2 -is the distance between the transfered proton and the axial O -CH 2 -group. O-H-O angle is the angle of the proton transfer, in degrees.H -O X H -O P H -O -CH 2 - O-H-O H -O X H -O P H -O -CH 2 - O-H-O

			HOH + EPH			CH3OH + EPH
	T S1 dd	1.308	1.147		133.3	1.255	1.185	131.7
	T S1 ud	1.199	1.269		129.2	1.196	1.265	127.6
	T S2 du		1.202	1.247	133.1		1.177	1.273	135.2
	T S2 uu		1.112	1.376	138.7		1.093	1.413	138.8
			HOH + REPH			CH3OH + REPH
	T S1 Exo du	1.199	1.272		129.0	1.197	1.266	127.5
	T S1 Exo dd	1.201	1.148		132.5	1.259	1.182	131.8
	T S1 Endo dd	1.298	1.152		132.7	1.255	1.185	131.7
	T S1 Endo ud	1.199	1.271		129.0	1.196	1.265	127.5
	T S2 Exo ud		1.205	1.243	133.9		1.205	1.243	133.9
	T S2 Exo uu		1.118	1.367	140.1		1.094	1.413	138.6
	T S2 Endo du		1.218	1.232	133.2		1.183	1.265	135.9
	T S2 Endo uu		1.118	1.367	140.1		1.099	1.402	140.4

5 Acknowledgment DY is grateful for financial support provided by the National Institutes of Health (Grant 1R01-GM62248-01A1). A.D. is grateful for financial support provided by the Centre National de la Recherche Scientifique (CNRS), the Institut National de la Sante et de la Recherche Medicale (INSERM), and the Universite L. Pasteur de Strasbourg. Computational resources were provided in part by the Minnesota Supercomputing Institute and in part by the Centre Informatique National de l'Enseignement Superieur (CINES) and the Institut du Developpement et des Ressources en Informatique Scientifique (IDRIS). The part of the work done at Harvard University was supported by the National Institutes of Health.

Associativeness Index

Figure 9: Relative aqueous free energies with respect to separate reactants as a function of a reaction coordinate (R.C.) defined as P -OX P -OX P -P -OC P -OC R . P-OC and P-OX are the two axial P-O interatomic distances, P-OC is the P-O distance with the leaving group and P-OX with the nucleophile. The two reference values P -OC R and P -OX P are the corresponding phosphate bonds in the reactant (R) and product (P), respectively. This reaction coordinate can be interpret as the difference between the the degree of phosphorus-nucleophile bond formation and bond cleavage with the leaving group for a given structure. Bottom diagram, represents the degree of associativeness of the structure as defined as 1 2 * P -OX P -OX P + P -OC P -OC R . A value of one represent the maximum degree of associativeness, and values larger than one indicate a higher dissociative character of the structure.