Modelling Complex Bimolecular Reactions in a Condensed Phase: The Case of Phosphodiester Hydrolysis

Abstract

(1) Background: the theoretical modelling of reactions occurring in liquid phase is a research line of primary importance both in theoretical–computational chemistry and in the context of organic and biological chemistry. Here we present the modelling of the kinetics of the hydroxide-promoted hydrolysis of phosphoric diesters. (2) Method: the theoretical–computational procedure involves a hybrid quantum/classical approach based on the perturbed matrix method (PMM) in conjunction with molecular mechanics. (3) Results: the presented study reproduces the experimental data both in the rate constants and in the mechanistic aspects (C–O bond vs. O–P bond reactivity). The study suggests that the basic hydrolysis of phosphodiesters occurs through a concerted $A_N{D}_N$ mechanism, with no formation of penta-coordinated species as reaction intermediates. (4) Conclusions: the presented approach, despite the approximations, is potentially applicable to a large number of bimolecular transformations in solution and therefore leads the way to a fast and general method to predict the rate constants and reactivities/selectivities in complex environments.

1. Introduction

Phosphodiester bonds are the linkages that connect biological information stored in nucleic acids. Their resistance to hydrolysis [1] makes them ideal for preserving biochemical information. However, phosphodiester hydrolysis is a crucial process in biochemical systems and for this reason, nature has developed a number of enzymes able to catalyse the hydrolysis of these compounds in aqueous solution. The rate enhancements induced by these enzymes (e.g., kinases, ATPases, phosphatases) are astonishing (up to ${10}^{17}$ in the case of staphylococcal nuclease) [2,3]. On the other hand, numerous research groups have designed and synthesized small molecule enzyme mimics able to cleave DNA, RNA and their model compounds in an effective and selective fashion [4,5,6,7].

The mechanism of phosphodiester hydrolysis was investigated at different pH conditions and for diverse phosphodiesters through experimental and theoretical–computational approaches for both catalysed [2,4,8,9,10,11,12,13,14] and uncatalysed reactions [9,15,16,17,18]. Both approaches present some critical aspects. The measurements of the rate constants for the uncatalysed hydrolysis of phosphodiesters are hampered by the extremely long half-life time of the reaction (e.g., 30,000,000 years at 25 °C) [2]. For this reason, there is the need to carry out kinetic experiments at elevated temperatures and pressures and obtain the reaction rate at room temperature by extrapolation [2]. Regarding the computational approach, the problems manifest in the choice of the proper method to describe the reaction coordinate and the solvent effect. For the reasons illustrated above, the development of an efficient and reliable method to investigate phosphodiester hydrolysis is crucial to provide useful insights into the mechanisms of these reactions and to evaluate the catalytic enhancement offered by the enzymes and artificial catalysts with respect to spontaneous hydrolysis.

In this work, we focused on the study of the hydrolysis of a simple phosphodiester, dineopentyl phosphate, hereafter referred to as $N{p}_2{P}^{-}$. The hydrolysis of this substrate was investigated by Wolfenden et al. through an in-depth kinetic analysis [2], and therefore experimental data are available for a direct comparison. This compound has the tendency, similarly to DNA, to undergo cleavage of the phosphorous–oxygen bond, instead of the carbon–oxygen bond as many phosphoric diesters [2]. This reactive behaviour can be attributed to the steric hindrance provided by the t-butyl units that force the attack of the upcoming nucleophile to the phosphorous atom, as illustrated and discussed in detail in the present manuscript (vide infra. Experimental data indicate that the hydroxide-catalyzed hydrolysis dominates over the spontaneous reaction as the leaving groups become poorer [2] and, according to this evidence, we have selected the hydroxide ion attack onto the substrate, instead of a water molecule, as a test reaction. The hydroxide attack is also frequent in enzyme-catalysed reactions, where a metal cation helps the deprotonation of a water molecule [8,19,20]. This transformation was modelled by a statistical-mechanical approach based on a quantum mechanics/molecular mechanics hybrid method. The results obtained by this approach are illustrated and discussed in detail and compared with the experimental data.

2. Theory

2.1. General Considerations

Let us consider a simple reaction from the reactant (R) chemical state to the product (P) chemical state, without involving any distinguishable intermediate chemical state, as occurring within a single reactive centre (RC), either a molecule or a complex, interacting with its environment. Note that we assume for the R and P chemical states a fully in equilibrium: i.e., the corresponding mean residence time is large enough to allow for internal relaxation. We can also reasonably assume that the reaction events must pass trough a specific and unique transition state ($TS$) with a mean traversing time too short for any internal relaxation, hence providing the reaction scheme

(1)$R\stackrel{k_R}{\to }TS\stackrel{k}{\to }P$

and thus

(2)$\begin{array}{ccc}\hfill \dot{\left[R\right]}& =& -k_R\left[R\right]\hfill \end{array}$

(3)$\begin{array}{ccc}\hfill \dot{\left[TS\right]}& =& k_R\left[R\right]-k\left[TS\right]\hfill \end{array}$

(4)$\begin{array}{ccc}\hfill \dot{\left[P\right]}& =& k\left[TS\right]\hfill \end{array}$

Considering that we necessarily have $k\gg k_R$ we can assume the steady state for $TS$ (i.e., $\dot{\left[TS\right]}\cong 0$) providing

(5)$\begin{array}{ccc}\hfill \left[TS\right]& \cong & \frac{k_R}{k}\left[R\right]\hfill \end{array}$

(6)$\begin{array}{ccc}\hfill \dot{\left[P\right]}& \cong & k\frac{k_R}{k}\left[R\right]=k_R\left[R\right]\hfill \end{array}$

It is reasonable to assume that within the $TS$ equilibrium ensemble the outward flux rate constants are identical, given by $k_{TS}\cong k/2$ (only one half of the $TS$ population is moving towards P), and thus from Equation (6) we have

(7)$\begin{array}{ccc}\hfill \dot{\left[P\right]}& \cong & -\dot{\left[R\right]}\cong \frac{k_R}{k_{TS}}\frac{k}{2}\left[R\right]=\frac{Q_{TS}}{Q_R}\frac{k}{2}\left[R\right]\hfill \\ & =& {\mathcal{K}}_R\left[R\right]\hfill \end{array}$

with $Q_{TS}$ and $Q_R$ the canonical partition functions of the whole system with a selected single RC complex (the reference RC) constrained in the $TS$ or the R state, respectively (note that for reactive centres at infinite dilution the whole system can be considered as including a single RC embedded into a huge amount of solvent molecules). From Equation (7) it readily follows that the reaction process can be described by the simplified reaction scheme

(8)$R\stackrel{{\mathcal{K}}_R}{\to }P$

2.2. The Reaction Coordinate, the Landau Free Energy and the Rate Constant

Let us assume that we can identify a generalized coordinate $\xi$ (i.e., the reaction coordinate) univocally defining, via its different intervals, the R, $TS$ and P chemical states of the reference RC (infinite dilution at the reactive centre the single RC is present), and thus properly describing the reaction process (see Figure 1). Defining with $[{\xi }_{TS}-\delta /2,{\xi }_{TS}+\delta /2]$ the $TS$ tiny reaction coordinate range and with ${\xi }_m$ the minimum $\xi$ value of the R chemical state, we can express the partition functions $Q_R$ and $Q_{TS}$ via the Landau free energy $\mathcal{A}$ defined by

(9)$\mathcal{A}\left(\xi \right)=-k_{B}Tln\mathcal{Q}\left(\xi \right)$

with $\mathcal{Q}\left(\xi \right)$ the partition function density providing

(10)$\begin{array}{ccc}\hfill Q_R& =& {\int }_{{\xi }_m}^{{\xi }_{TS}-\delta /2}\mathcal{Q}\left(\xi \right)d\xi ={\int }_{{\xi }_m}^{{\xi }_{TS}-\delta /2}{e}^{-\beta \mathcal{A}\left(\xi \right)}d\xi \cong {\int }_{{\xi }_m}^{{\xi }_{TS}}{e}^{-\beta \mathcal{A}\left(\xi \right)}d\xi \hfill \end{array}$

(11)$\begin{array}{ccc}\hfill Q_{TS}& =& {\int }_{{\xi }_{TS}-\delta /2}^{{\xi }_{TS}+\delta /2}\mathcal{Q}\left(\xi \right)d\xi ={\int }_{{\xi }_{TS}-\delta /2}^{{\xi }_{TS}+\delta /2}{e}^{-\beta \mathcal{A}\left(\xi \right)}d\xi \cong e^{-\beta \mathcal{A}\left({\xi }_{TS}\right)}\delta \hfill \end{array}$

where $1/\beta =k_{B}T$ and we used $\mathcal{A}\left(\xi \right)\cong \mathcal{A}\left({\xi }_{TS}\right)$ within the $TS$ range.

From Equations (7), (10) and (11) we obtain

(12)$\begin{array}{ccc}\hfill {\mathcal{K}}_R& \cong & \frac{e^{-\beta \mathsf{\Delta }{\mathcal{A}}_R^{†}}\delta }{{\int }_{{\xi }_m}^{{\xi }_{TS}}{e}^{-\beta \mathsf{\Delta }\mathcal{A}\left(\xi \right)}d\xi }\frac{k}{2}\hfill \end{array}$

(13)$\begin{array}{ccc}\hfill \mathsf{\Delta }\mathcal{A}\left(\xi \right)& =& \mathcal{A}\left(\xi \right)-\mathcal{A}\left({\xi }_R\right)\hfill \end{array}$

(14)$\begin{array}{ccc}\hfill \mathsf{\Delta }{\mathcal{A}}_R^{†}& =& \mathcal{A}\left({\xi }_{TS}\right)-\mathcal{A}\left({\xi }_R\right)\hfill \end{array}$

where ${\xi }_R$ is a reference position within the R state range (typically the Landau free energy minimum). Realizing that the $TS$, corresponding to a tiny $\xi$ interval, is characterized by a constant density within the whole $\delta$ range and for each reaction event by a virtually constant $d\xi /dt$ velocity, we can express the outward flux (forming the P population) via $k\left[TS\right]\cong \left[TS\right]⟨vs.⟩/\delta$ with $⟨vs.⟩$ the mean $TS$ traversing velocity given by the equilibrium mean value of $v=|d\xi /dt|$. Therefore, it follows

(15)$k\cong \frac{⟨vs.⟩}{\delta }$

readily providing

(16)${\mathcal{K}}_R\cong \frac{e^{-\beta \mathsf{\Delta }{\mathcal{A}}_R^{†}}}{{\int }_{{\xi }_m}^{{\xi }_{TS}}{e}^{-\beta \mathsf{\Delta }\mathcal{A}\left(\xi \right)}d\xi }\frac{⟨vs.⟩}{2}$

Finally, we can obtain the Landau free energy change along $\xi$, via [21]

(17)$\begin{array}{ccc}\hfill \mathcal{A}\left(\xi \right)-\mathcal{A}\left({\xi }_{ref}\right)& =& A_{\xi }-A_{{\xi }_{ref}}\cong -k_{B}Tln{⟨e^{-\beta \mathsf{\Delta }\mathcal{U}}⟩}_{{\xi }_{ref}}\hfill \end{array}$

(18)$\begin{array}{ccc}\hfill A_{\xi }& =& \mathcal{A}\left(\xi \right)-k_{B}Tln\delta \hfill \end{array}$

(19)$\begin{array}{ccc}\hfill A_{{\xi }_{ref}}& =& \mathcal{A}\left({\xi }_{ref}\right)-k_{B}Tln\delta \hfill \end{array}$

(20)$\begin{array}{ccc}\hfill \mathsf{\Delta }\mathcal{U}& =& \mathcal{U}\left(\xi \right)-\mathcal{U}\left({\xi }_{ref}\right)\hfill \end{array}$

where the subscript of the angle brackets indicate averaging within the ${\xi }_{ref}$ reference ensemble, $\mathcal{U}\left(\xi \right)$ and $\mathcal{U}\left({\xi }_{ref}\right)$ are the perturbed electronic energies (estimated by means of the perturbed matrix method (PMM) [22,23]) as obtained either at fixed $\xi$ or at fixed ${\xi }_{ref}$ for each accessible semi-classical configuration, once relaxed all the quantum vibrational coordinates are at their minimum energy position (we consider the quantum vibrational energy as independent of the $\xi$ coordinate), and $A_{\xi }$ and $A_{{\xi }_{ref}}$ are the Helmholtz free energies when confining the RC within the tiny $\delta$ range either around $\xi$ or around ${\xi }_{ref}$. Note that in practice Equation (17) is accurate only when $\xi$ is not too far from ${\xi }_{ref}$, thus allowing us to use the fixed ${\xi }_{ref}$ MD simulation to reconstruct the $\xi$ ensemble (i.e., the MD simulation provides a reasonable sampling of the relevant phase space regions of both ensembles). Therefore, for evaluating the complete free energy profile we can use a set of reference positions ${\xi }_{ref,n},n=1,N$ providing a set of $[{\xi }_{ref,n},{\xi }_{ref;n+1}]$ sub-ranges each to be used to obtain the corresponding (inner) Landau free energy changes. By summing the free energy variations over such sub-ranges the complete Landau free energy profile can be reconstructed and hence by averaging over the free energy profiles as obtained by both the ${\xi }_{ref,1}\to {\xi }_{ref,N}$ and ${\xi }_{ref,N}\to {\xi }_{ref,1}$ calculations we can obtain a reliable estimate of the reaction Landau free energy surface.

2.3. Application to the $N{p}_{2}P$ Hydrolysis Reaction

The general scheme for the hydrolysis reaction considered is

(21)$N{p}_2{P}^{-}+H{O}^{-}\underset{k_2}{\overset{k_A}{\underset{\leftharpoondown }{\rightharpoonup }}}{C}_{NR}\underset{k_1}{\overset{k_{-1}}{\underset{\leftharpoondown }{\rightharpoonup }}}{C}_R\stackrel{{\mathcal{K}}_R}{\to }P$

where $N{p}_2{P}^{-}$ is the dineopentyl phosphate, $H{O}^{-}$ is the hydroxide ion and $C_{NR},C_R$ are the non-reactive and reactive complexes, respectively: the former defined by the interacting partners still insufficiently close to alter their covalent structure while the latter corresponding to closely interacting partners possibly undergoing covalent rearrangement (i.e., the R state of Scheme (8)). In the previous reaction scheme the association rate constant $k_A$ can be obtained via [24,25]

(22)$k_A=4\pi (D_{N{p}_2{P}^{-}}+D_{H{O}^{-}})r_0{\mathcal{N}}_A$

where $r_0$ is the maximum distance between $N{p}_2{P}^{-}$ and $H{O}^{-}$ within the non-reactive complex (i.e., the radius length of the sphere defining the overall complex $C=C_{NR}+C_R$), $D_{N{p}_2{P}^{-}}$ and $D_{H{O}^{-}}$ are the free diffusion coefficients of the two reactant molecules characterizing their behaviour for $r\ge r_0$, and ${\mathcal{N}}_A$ is the Avogadro number. Finally, assuming steady state for $C_{NR}$ and $C_R$, we have [25]

(23)$\dot{\left[N{p}_2{P}^{-}\right]}\cong -k_H\left[N{p}_2{P}^{-}\right]\left[H{O}^{-}\right]$

with

(24)$k_H=\frac{{\mathcal{K}}_R{k}_A{k}_{-1}}{(k_2+k_{-1})(k_1+{\mathcal{K}}_R)-k_1{k}_{-1}}\cong \frac{k_A}{k_2}\frac{k_{-1}}{k_1}{\mathcal{K}}_R$

as the overall (hydrolysis) rate constant (note that in Equation (24) we used $k_1+{\mathcal{K}}_R\cong k_1$ as following from ${\mathcal{K}}_R\ll k_1$, see Table 1). By means of proper MD simulations, the PMM approach and the relations described in the theory section we obtained all the kinetic constants to be used in Equation (24).

3. Methods

3.1. Quantum Chemical Calculations

The quantum chemical calculations were conducted using Orca software package [26]. To model the hydrolysis of dineopentyl phosphate ($N{p}_2{P}^{-}$) with the attack on the phosphorous atom, we selected as the quantum centre (QC) the dimethyl-substituted phosphate ($DM{P}^{-}$) and the hydroxide ion. The quantum chemical calculations were performed using density functional theory (DFT) to obtain the intrinsic reaction coordinate (IRC). The IRC calculation requires a guess of the transition state (TS) of the reaction. The initial guess for the transition state structure was taken from previous work [8] and then optimized for a saddle point at the DFT/PBE0 [27] level of theory using the 6-311++G(2d,2p) [28,29] basis set. The analysis of the transition state hessian confirmed that the obtained structure corresponds to a saddle point of the potential energy surface (PES).

The choice of the PBE0 functional was motivated by a previous benchmark work present in the literature [30].

To show that the attack of the hydroxide ion on the neopentyl group is energetically precluded, we performed a PES scan along the hydroxide oxygen atom–dineopentyl carbon atom and the dineopentyl carbon atom–leaving group (see Figure 2a) distances, at DFT/PBE0 [27] level of theory using the 6-31+G(d,p) [28,29,31] basis set. The Cartesian coordinates of the transition states found and used in the present work are available in the Supplementary Materials.

3.2. Molecular Dynamics Simulations

The molecular dynamics (MD) simulations were conducted starting from the structures obtained from the IRC calculation (as explained in the previous subsection) and by replacing one hydrogen atom of each methyl group of the $DM{P}^{-}$ with a t-butyl group to obtain a model of the dineopentyl phosphate and hydroxide ion, along the reaction coordinate.

All the molecular dynamics simulations have been performed using Gromacs 2022 software package [32]. The initial structure was centred in a cubic box and solvated using the SPC [33] water model, whereas the force field parameters of the neopentyl phosphate were obtained by means of the ACPYPE package [34]. Two ${\mathrm{Na}}^{+}$ ions have been inserted in the simulation box to reach the electro-neutrality. An energy minimization step was performed using the steepest descent algorithm freezing the atoms corresponding to the QC. After the minimization, a series of equilibration steps, lasting 50 ps with a timestep of 2 fs, were performed in the NVT ensemble: after each run, we checked the density of the system and tune the size of the box in order to reach the correct value of density. After the equilibration steps, each system has been simulated for 25 ns at fixed temperature and volume. The temperature was kept constant at 300 K using the V-rescale thermostat [35] using a ${\tau }_T$ of 0.1 ps. The electrostatic interactions were calculated using the particle mesh Ewald method [36,37] with a cut-off of 1.2 nm. A cut-off of 1.2 nm was used for the van der Waals interactions.

Because the reaction shows a relevant reorganization of the internal degrees of freedom, we divided the intrinsic reaction coordinate into five subparts. For each of them, an MD simulation (maintaining the QC fixed, in a reference structure for the selected subpart) was carried out to provide a proper semi-classical perturbation of the QC.

Additionally, we performed an MD simulation, lasting 100 ns, of the free reactants using a timestep of 2 fs to obtain the hydroxide ion–dineopentyl phosphate distance probability distribution.

A short simulation, lasting 200 ps, with a timestep of 1 fs of the dineopentyl phosphate in aqueous solution was performed to estimate its diffusion coefficient according to the Einstein Equation [38]. For the hydroxide ion, we used the experimental diffusion coefficient reported in [39]. The diffusion coefficients are required for the calculation of the kinetic constant $k_A$ (see Scheme (21)) according to Equation (22).

We also performed a set of 70 simulations, lasting 500 ps, starting from random initial configuration within the hydroxide ion–dineopentyl phosphate distance of 0.5–0.6 nm to reconstruct the reactive and non-reactive complex dissociation kinetics.

3.3. Computational Strategy

The computational strategy can be summarized in the following steps:

• Calculation of the unperturbed quantum chemical properties, needed to apply the MD-PMM procedure, for each point along the reaction coordinate (electronic ground state energy, permanent dipole moments and CHELPG charges [40]).

• Simulations of selected representative ensemble along the reaction coordinate (reactants, two intermediate structures between reactants and transition state, corresponding to different nucleophile–phosphorous atom distances, and the transition state) providing the perturbation along the reaction coordinate.

• MD-PMM calculation for each point to obtain the free energy difference from one point to the following along the reaction coordinate, according to Equation (17).

• Calculation of the kinetic rate constants involved in the reactive non-reactive equilibrium in Scheme (21) by performing different sets of MD simulations: (i) a 100 ns long simulation of the free reactants to obtain the hydroxide ion–dineopentyl phosphate distance probability distribution and the equilibrium fraction of the non-reactive complex, ${\chi }_{eq}$; (ii) one short simulation for the determination of the diffusion coefficient of the dineopentyl phosphate; (iii) a set of 70 simulations in order to reconstruct the reactive and non-reactive complex dissociation kinetics.

4. Results and Discussion

As explained in the Theory section, the evaluation of the overall rate constant $k_H$ for the P–O cleavage reaction of the $N{p}_2{P}^{-}$ in water requires the estimate of the rate constants associated with the three reaction steps, described in Scheme (21), using Equation (24).

4.1. Evaluation of the Rate Constant ${\mathcal{K}}_R$ of the Dineopentyl Phosphate Hydrolysis

The methodology described in the Theory and Methods sections was applied to reconstruct the full kinetics of the hydrolysis of the dineopentyl phosphate. This substrate was selected because of the availability of experimental data [2], as highlighted in the introduction. It is known that the hydrolysis of the dialkyl phosphodiesters with minor steric hindrance around the carbon atoms bound to the oxygen atoms proceeds preferentially via the nucleophilic attack at this carbon atom with respect to the phosphorous atom [17,18] (Figure 2a). However, in the case of the dineopentyl phosphate, the steric hindrance drives the attack onto the phosphorous atom [1], as shown in Figure 2b.

Therefore, to quantitatively model the free energy profile of the reaction by means of the PMM-MD approach, we first evaluate the minimum energy path along both the P–O and C–O bond cleavages. To study the attack of the hydroxide ion on the phosphorous atom of the dineopentyl phosphate, we performed a TS optimization of the quantum centre and analysis of the TS hessian to confirm that the obtained structure is a saddle point of the potential energy surface of the system. From this TS structure, an intrinsic reaction coordinate calculation was performed to obtain the reaction profile reported in Figure 3b. In the case of C–O bond cleavage, the same calculation on $N{p}_2{P}^{-}$ failed because of the large steric hindrance due to the neopentyl groups. Therefore, a bi-dimensional PES along the two bond distances involved in the reaction, i.e., the distance between the carbon atom and the oxygen of the hydroxide ($d_{C-OH}$) and between that carbon atom and the oxygen of the leaving group ($d_{C-OLg}$), was performed. These results reported in Figure 3 clearly show a higher energy barrier associated with the C–O bond cleavage (36.1 kcal/mol) with respect to the P–O bond (27.5 kcal/mol). This indicates that the steric hindrance inhibits the C–O bond cleavage. Therefore, we have applied the MD-PMM procedure to evaluate only the thermodynamics and kinetics in solution for the P–O bond cleavage.

To this end, the electronic energies were evaluated for ∼100 structures taken along the reaction coordinate and perturbed—in the corresponding MD ensemble—by the electric field as provided by the MD simulations. These data indicate that the effect of the environment is to stabilize both the reactants and the TS, as shown by the distributions of the reactant and TS perturbed electronic energies (see Figure 4).

However, the TS perturbed energies decrease less than the reactant state ones, resulting in a higher barrier with respect to the vacuum condition (see Table 2). This effect can be ascribed to a better ability of polar solvent molecules to stabilize two negative charges located on two separated molecular entities compared with a doubly charged molecule.

Using these perturbed energy trajectories in Equation (17), it was possible to calculate the corresponding (Helmoltz) free energy profile along the reaction coordinate, which results in a barrier of 38.8 kcal/mol, Figure 5. Using the statistical–mechanical approach described in the Theory section, i.e., Equation (16), we obtain, for the last (and slowest) reaction step (see Scheme (21)), a kinetic rate constant ${\mathcal{K}}_R$ as high as 1.85·${10}^{-15}$ s${}^{-1}$.

4.2. Evaluation of the Other Kinetic Rate Constants

In order to evaluate the other rate constants for the reaction steps providing the complex dissociation and the non-reactive reactive complex transitions, we computed the overall complex dissociation rate constant ($k_D={\chi }_{eq}{k}_2,{\chi }_{eq}=k_1/(k_1+k_{-1})$), assuming pre-equilibrium within the reactant complex, as provided by 70 MD simulations starting from random initial configurations representative of the reactive complex. The direct evaluation of $k_D$ can be obtained by monitoring and fitting the overall complex fraction in time, i.e., by first-order modelling the time decay of the trajectory fraction corresponding to the reactants with distances within 1.2 nm. This cut-off value, defining the overall complex radius $r_0$, can be identified by inspecting the hydroxide ion–dineopentyl phosphate distance distribution (Figure 6) as provided by MD simulations. When that distribution starts to exhibit a purely quadratic behaviour, the two reactants can be considered to be outside of the non-reactive complex region.

Analogously, the rate constant $k_1$ for the reactive complex → non-reactive complex transition can be computed by monitoring the time decay of the complex fraction corresponding to the reactants with distances within 0.6 nm (i.e., the cut-off radius defining the reactive complex), see Figure 7. The computed values for $k_D$ and $k_1$ are reported in Table 1.

Moreover, from the 100 ns MD simulation of the free reactants we obtained the equilibrium fraction of the non-reactive complex (${\chi }_{eq}$) used to estimate the kinetic constants $k_{-1}$ and $k_2$ (see Scheme (21)), also reported in Table 1.

Finally, using all the rate constants for the cleavage reaction steps in Equation (24), it was possible to estimate the overall rate constant $k_H=2.61$·${10}^{-15}$ s${}^{-1}$ in very good agreement with the experimental data, both reported in Table 1.

Our data clearly show that the first two steps of the reaction have a limited effect on the overall rate of the reaction, as pointed out by the limited difference between ${\mathcal{K}}_R$ and $k_H$.

4.3. Hydrolysis of the Dimethyl Phosphate

In order to confirm that the steric hindrance is able to drive the reaction towards different paths, that is P–O or C–O bond cleavage, the slow step of the hydrolysis reaction of the dimethyl phosphate ($DM{P}^{-}$) was investigated by means of the same theoretical–computational strategy for both C–O and P–O bond breaks. In agreement with the experimental data, the free energy barrier is remarkably lower for the C–O bond break with respect to the cleavage of the P–O bond.

Assuming that in the case of the $DM{P}^{-}$ the first two reaction steps exhibit similar behaviours with respect to $N{p}_2{P}^{-}$, we compare the ratio between the experimental overall rate constants of the two possible bond breaks with the ratio of the corresponding ${\mathcal{K}}_R$. From these data, it comes out that our procedure correctly predicts the favourite path (see Table 2).

4.4. Mechanistic Considerations

The hydrolysis of phosphodiesters may proceed via a step-wise mechanism or, alternatively, a concerted pathway (see Figure 8). In the first case, the formation of the intermediate can be generated either by the attack of the nucleophile onto the phosphate unit, i.e., $A_N+D_N$, mechanism with the generation of a penta-coordinated species, or by the early departure of the leaving group (Lg), named as $D_N+A_N$ [8,9,41,42]. The latter case can be considered only in the case of substrates provided with very good leaving groups. On the other hand, in a concerted reaction mechanism, the penta-coordinated phosphorane is a transition state instead of an intermediate, namely $A_N{D}_N$ pathway, indicated also as $S_{N}2$(P) for its close similarity to carbon nucleophilic substitution. Even in the concerted transformation, the TS can verge to one of the two extreme conditions illustrated above. A transition state that is closer to the associative limit is labelled as $A_N{D}_N$, i.e., asynchronous tight. On the contrary, an asynchronous loose $D_N{A}_N$ mechanism is closer to the dissociative limit [8,41]. From the analysis of the potential energy surface and the corresponding O${}_{nucleophile}$–P and P–OLg distances (see Figure 3) the reaction pathway does not indicate the existence of any intermediate, as no energy minima is detected along the profile, and therefore the mechanism can be considered concerted and synchronous, or slightly “associative” ($A_N{D}_N$). This is in contrast to what has been observed in a previous investigation based on a more approximate theoretical–computational approach and on the implicit treatment of the solvent [9].

5. Conclusions

Our theoretical and computational procedure, based on a mixed quantum/classical approach, was employed to investigate a chemical transformation of relevant importance in biological systems. The proposed methodology, based on the perturbed matrix method (PMM), has turned out to be effective in the quantitative prediction of the kinetics and thermodynamics of the hydrolysis of phosphate diesters in water solution. The results obtained by such an approach, in spite of the adopted approximation, are in excellent agreement with the experimental rate constants. Furthermore, this approach provides useful mechanistic insights about the reaction pathway, i.e., C–O bond cleavage vs. O–P bond cleavage for phosphodiesters featuring different steric hindrance conditions, in accordance with the experimental picture. The analysis of the potential energy surface does not show the presence of any energy minimum along the reaction pathway, thus suggesting that the hydroxide-promoted hydrolysis of these compounds proceeds through a concerted “synchronous” mechanism. We believe the present approach can be applied in a wide variety of chemical transformations and therefore can open up the way to a fast and reliable general method to treat bimolecular reactions in condensed phase and in other complex environments.

Author Contributions

Writing-original draft: A.N.N., R.S., Writing-review and editing: R.S., M.D., Software, Writing-original draft: A.N.N., Formal analysis: A.N.N., A.O., Conceptualization: A.A., M.D., Methodology: A.A., Project administration: M.D. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Pictorial representation of a generic free energy profile along the reaction coordinate [Forumla omitted. See PDF.].

Figure 2. Schematic representation of the alternative reaction mechanisms for the cleavage of [Forumla omitted. See PDF.]: (a) C–O cleavage, [Forumla omitted. See PDF.]; and (b) P–O cleavage, [Forumla omitted. See PDF.](P) or [Forumla omitted. See PDF.].

Figure 3. (a) Bi-dimensional scan performed along the hydroxide oxygen atom–carbon atom and the carbon atom –leaving group oxygen distances; (b) Gas phase electronic energy profile along the reaction coordinate [Forumla omitted. See PDF.] (defined as the difference between the HO[Forumla omitted. See PDF.]–P and the P–OLg distances).

Figure 4. Reactants and TS perturbed electronic energies distributions with respect to their corresponding unperturbed states.

Figure 5. Free energy profile along the reaction coordinate [Forumla omitted. See PDF.] defined as the difference between the HO[Forumla omitted. See PDF.]–P and the P–OLg distances, for the reaction in aqueous solution.

Figure 6. Distribution of hydroxide ion–dineopentyl phosphate distance (navy); quadratic fitting function of the distance distribution (orange).

Figure 7. Time decay of the (total) complex (a) and the reactive complex (b) in blue dots and their corresponding exponential fits in solid black lines.

Figure 8. More O’Ferrall–Jencks diagram for the hydrolysis of phosphodiesters. Definition of the associative and dissociative mechanisms.

Table 1

Kinetic rate constants of the P–O bond cleavage in $N{p}_2{P}^{-}$.

Calculated Rate Constants	Value	Dimension
$k_A$	$0.056\pm 0.009$	${\mathrm{ps}}^{-1}{\mathrm{M}}^{-1}$
$k_D$	$0.0056\pm 0.0009$	${\mathrm{ps}}^{-1}$
$k_1$	$0.030\pm 0.002$	${\mathrm{ps}}^{-1}$
$k_{-1}$	$0.0049\pm 0.0003$	${\mathrm{ps}}^{-1}$
$k_2$	$0.007\pm 0.001$	${\mathrm{ps}}^{-1}$
${\mathcal{K}}_R$	$(1.85\pm 0.92)·{10}^{-15}$	${\mathrm{s}}^{-1}$
$k_H$	$(2.61\pm 1.23)·{10}^{-15}$	${\mathrm{s}}^{-1}{\mathrm{M}}^{-1}$
$k_H{\left[H{O}^{-}\right]}^{\ominus }$	$(2.61\pm 1.23)·{10}^{-15}$	${\mathrm{s}}^{-1}$
Experimental rate constant
$k_H$	1.4 $·{10}^{-15}$	${\mathrm{s}}^{-1}{\mathrm{M}}^{-1}$
$k_H{\left[H{O}^{-}\right]}^{\ominus }$ ${}^a$	1.4 $·{10}^{-15}$	${\mathrm{s}}^{-1}$

^a Hydroxide ion concentration ${\left[H{O}^{-}\right]}^{\ominus }$ = 1 M.

Table 2

Comparison of the hydrolysis of $DM{P}^{-}$ and $N{p}_2{P}^{-}$. $\mathsf{\Delta }{\mathcal{U}}^0$ corresponds to the unperturbed electronic energy barrier and $\mathsf{\Delta }A$ to the free energy barrier.

	${\mathit{DMP}}^{-}$		${\mathit{Np}}_2{\mathit{P}}^{-}$
Bond Cleavage	P–O	C–O	P–O	C–O
$\mathsf{\Delta }{\mathcal{U}}^0$ (kcal/mol)	27.5	28.0	28.5	36.1
$\mathsf{\Delta }A$ (kcal/mol)	39.5	30.9	38.8	-
${\mathcal{K}}_R$ (s${}^{-1}$)	3.75 $·{10}^{-16}$	3.17 $·{10}^{-10}$	1.85 $·{10}^{-15}$	-
$k_H^{exp}$ (s${}^{-1}$M${}^{-1}$)	NA	3.0 $·{10}^{-11}$	1.4 $·{10}^{-15}$	NA

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/molecules28052152/s1.

References

1. Wolfenden, R. Degrees of difficulty of water-consuming reactions in the absence of enzymes. Chem. Rev.; 2006; 106, pp. 3379-3396. [DOI: https://dx.doi.org/10.1021/cr050311y]

2. Schroeder, G.K.; Lad, C.; Wyman, P.; Williams, N.H.; Wolfenden, R. The time required for water attack at the phosphorus atom of simple phosphodiesters and of DNA. Proc. Natl. Acad. Sci. USA; 2006; 103, pp. 4052-4055. [DOI: https://dx.doi.org/10.1073/pnas.0510879103]

3. Zalatan, J.G.; Herschlag, D. Alkaline phosphatase mono-and diesterase reactions: Comparative transition state analysis. J. Am. Chem. Soc.; 2006; 128, pp. 1293-1303. [DOI: https://dx.doi.org/10.1021/ja056528r]

4. Raynal, M.; Ballester, P.; Vidal-Ferran, A.; van Leeuwen, P.W. Supramolecular catalysis. Part 2: Artificial enzyme mimics. Chem. Soc. Rev.; 2014; 43, pp. 1734-1787. [DOI: https://dx.doi.org/10.1039/C3CS60037H]

5. Mikkola, S.; Lönnberg, T.; Lönnberg, H. Phosphodiester models for cleavage of nucleic acids. Beilstein J. Org. Chem.; 2018; 14, pp. 803-837. [DOI: https://dx.doi.org/10.3762/bjoc.14.68]

6. Mancin, F.; Scrimin, P.; Tecilla, P. Progress in artificial metallonucleases. Chem. Commun.; 2012; 48, pp. 5545-5559. [DOI: https://dx.doi.org/10.1039/c2cc30952a]

7. Morrow, J.R.; Amyest, T.L.; Richard, J.P. Phosphate Binding Energy and Catalysis by Small and Large Molecules. Acc. Chem. Res.; 2008; 41, pp. 539-548. [DOI: https://dx.doi.org/10.1021/ar7002013]

8. Segreto, G.E.; Alba, J.; Salvio, R.; D’Abramo, M. DNA cleavage by endonuclease I-DmoI: A QM/MM study and comparison with experimental data provide indications on the environmental effects. Theor. Chem. Acc.; 2020; 139, pp. 1-7. [DOI: https://dx.doi.org/10.1007/s00214-020-2585-0]

9. Salvio, R.; Casnati, A. Guanidinium Promoted Cleavage of Phosphoric Diesters: Kinetic Investigations and Calculations Provide Indications on the Operating Mechanism. J. Org. Chem.; 2017; 82, pp. 10461-10469. [DOI: https://dx.doi.org/10.1021/acs.joc.7b01925]

10. Salvio, R.; Moliterno, M.; Caramelli, D.; Pisciottani, L.; Antenucci, A.; D’Amico, M.; Bella, M. Kinetic resolution of phosphoric diester by Cinchona alkaloid derivatives provided with a guanidinium unit. Cat. Sci. Technol.; 2016; 6, pp. 2280-2288. [DOI: https://dx.doi.org/10.1039/C5CY01208B]

11. Salvio, R.; Casnati, A.; Mandolini, L.; Sansone, F.; Ungaro, R. ATP cleavage by cone tetraguanidinocalix[4]arene. Org. Biomol. Chem.; 2012; 10, pp. 8941-8943. [DOI: https://dx.doi.org/10.1039/c2ob26688a]

12. Salvio, R.; Volpi, S.; Cacciapaglia, R.; Sansone, F.; Mandolini, L.; Casnati, A. Phosphoryl Transfer Processes Promoted by a Trifunctional Calix[4]arene Inspired by DNA Topoisomerase I. J. Org. Chem.; 2016; 81, pp. 9012-9019. [DOI: https://dx.doi.org/10.1021/acs.joc.6b01643] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27579493]

13. Cacciapaglia, R.; Di Stefano, S.; Mandolini, L.; Salvio, R. Reactivity of carbonyl and phosphoryl groups at calixarenes. Supramol. Chem.; 2013; 25, pp. 537-554. [DOI: https://dx.doi.org/10.1080/10610278.2013.824578]

14. Salvio, R.; Volpi, S.; Folcarelli, T.; Casnati, A.; Cacciapaglia, R. A calix[4]arene with acylguanidine units as an efficient catalyst for phosphodiester bond cleavage in RNA and DNA model compounds. Org. Biomol. Chem.; 2019; 17, pp. 7482-7492. [DOI: https://dx.doi.org/10.1039/C9OB01141B] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31364652]

15. Imhof, P.; Fischer, S.; Krämer, R.; Smith, J.C. Density functional theory analysis of dimethylphosphate hydrolysis: Effect of solvation and nucleophile variation. J. Mol. Struct. Theochem; 2005; 713, pp. 1-5. [DOI: https://dx.doi.org/10.1016/j.theochem.2004.10.028]

16. Lopez, X.; Dejaegere, A.; Leclerc, F.; York, D.M.; Karplus, M. Nucleophilic attack on phosphate diesters: A density functional study of in- line reactivity in dianionic, monoanionic, and neutral systems. J. Phys. Chem. B; 2006; 110, pp. 11525-11539. [DOI: https://dx.doi.org/10.1021/jp0603942]

17. Kirby, A.J.; Younas, M. The Reactivity of Phosphate Esters. Diester Hydrolysis. J. Chem. Soc. B; 1970; pp. 510-513. [DOI: https://dx.doi.org/10.1039/j29700000510]

18. Kim, J.H.; Chin, J. Dimethyl Phosphate Hydrolysis at Neutral pH. J. Am. Chem. Soc.; 1992; 114, pp. 9792-9795. [DOI: https://dx.doi.org/10.1021/ja00051a011]

19. Prieto, J.; Epinat, J.C.; Redondo, P.; Ramos, E.; Padro, D.; Cedrone, F.; Montoya, G.; Paques, F.; Blanco, F. Generation and Analysis of Mesophilic Variants of the Thermostable Archaeal I-DmoI Homing Endonuclease. J. Biol. Chem.; 2008; 283, pp. 4364-4374. [DOI: https://dx.doi.org/10.1074/jbc.M706323200]

20. Stoddart, B.L. Homing endonuclease structure and function. Q. Rev. Biophys.; 2006; 38, pp. 49-95. [DOI: https://dx.doi.org/10.1017/S0033583505004063]

21. Amadei, A.; D’Alessandro, M.; Aschi, M. Statistical mechanical modeling of chemical reactions in complex systems: The reaction free energy surface. J. Phys. Chem. B; 2004; 108, pp. 16250-16254. [DOI: https://dx.doi.org/10.1021/jp048788l]

22. Zanetti-Polzi, L.; Del Galdo, S.; Daidone, I.; D’Abramo, M.; Barone, V.; Aschi, M.; Amadei, A. Extending the perturbed matrix method beyond the dipolar approximation: Comparison of different levels of theory. Phys. Chem. Chem. Phys.; 2018; 20, pp. 24369-24378. [DOI: https://dx.doi.org/10.1039/C8CP04190C] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30215645]

23. Chen, C.G.; Nardi, A.N.; Amadei, A.; D’Abramo, M. PyMM: An Open-Source Python Program for QM/MM Simulations Based on the Perturbed Matrix Method. J. Chem. Theory Comput.; 2023; 19, pp. 33-41. [DOI: https://dx.doi.org/10.1021/acs.jctc.2c00767] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36378163]

24. Cantor, C.R.; Schimmel, P.R. Biophysical Chemistry Part II; W.H. Freeman and Co.: San Francisco, CA, USA, New York, NY, USA, 1998.

25. Del Galdo, S.; Aschi, M.; Amadei, A. In silico characterization of bimolecular electron transfer reactions: The ferrocene–ferrocenium reaction as a test case. Int. J. Quantum Chem.; 2016; 116, pp. 1723-1730. [DOI: https://dx.doi.org/10.1002/qua.25212]

26. Neese, F. The ORCA program system. Wiley Interdiscip. Rev. Comput. Mol. Sci.; 2012; 2, pp. 73-78. [DOI: https://dx.doi.org/10.1002/wcms.81]

27. Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys.; 1999; 110, pp. 6158-6170. [DOI: https://dx.doi.org/10.1063/1.478522]

28. Krishnan, R.; Binkley, J.S.; Seeger, R.; Pople, J.A. Self-consistent molecular orbital methods. XX. A basis set for correlated wave functions. J. Chem. Phys.; 1980; 72, pp. 650-654. [DOI: https://dx.doi.org/10.1063/1.438955]

29. Clark, T.; Chandrasekhar, J.; Spitznagel, G.W.; Schleyer, P.V.R. Efficient diffuse function-augmented basis sets for anion calculations. III. The 3-21+ G basis set for first-row elements, Li–F. J. Comput. Chem.; 1983; 4, pp. 294-301. [DOI: https://dx.doi.org/10.1002/jcc.540040303]

30. Ribeiro, A.J.; Ramos, M.J.; Fernandes, P.A. Benchmarking of DFT functionals for the hydrolysis of phosphodiester bonds. J. Chem. Theory Comput.; 2010; 6, pp. 2281-2292. [DOI: https://dx.doi.org/10.1021/ct900649e]

31. Hehre, W.J.; Ditchfield, R.; Pople, J.A. Self—consistent molecular orbital methods. XII. Further extensions of Gaussian—type basis sets for use in molecular orbital studies of organic molecules. J. Chem. Phys.; 1972; 56, pp. 2257-2261. [DOI: https://dx.doi.org/10.1063/1.1677527]

32. Abraham, M.J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J.C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX; 2015; 1–2, pp. 19-25. [DOI: https://dx.doi.org/10.1016/j.softx.2015.06.001]

33. Mark, P.; Nilsson, L. Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K. J. Phys. Chem. A; 2001; 105, pp. 9954-9960. [DOI: https://dx.doi.org/10.1021/jp003020w]

34. Sousa da Silva, A.W.; Vranken, W.F. ACPYPE-Antechamber python parser interface. BMC Res. Notes; 2012; 5, pp. 1-8. [DOI: https://dx.doi.org/10.1186/1756-0500-5-367] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22824207]

35. Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys.; 2007; 126, 014101. [DOI: https://dx.doi.org/10.1063/1.2408420]

36. Darden, T.; York, D.; Pedersen, L.G. Particle mesh Ewald: An N-log(N) method for Ewald sums in large systems. J. Chem. Phys.; 1993; 98, pp. 10089-10092. [DOI: https://dx.doi.org/10.1063/1.464397]

37. Essmann, U.; Perera, L.; Berkowitz, M.L.; Darden, T.; Lee, H.; Pedersen, L.G. A smooth particle mesh Ewald method. J. Chem. Phys.; 1995; 103, pp. 8577-8593. [DOI: https://dx.doi.org/10.1063/1.470117]

38. Allen, M. DJ Tildesley Computer Simulation of Liquids; Clarendon: Oxford, UK, 1987.

39. Light, T.S.; Licht, S.; Bevilacqua, A.C.; Morash, K.R. The fundamental conductivity and resistivity of water. Electrochem. Solid-State Lett.; 2004; 8, E16. [DOI: https://dx.doi.org/10.1149/1.1836121]

40. Breneman, C.M.; Wiberg, K.B. Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis. J. Comput. Chem.; 1990; 11, pp. 361-373. [DOI: https://dx.doi.org/10.1002/jcc.540110311]

41. Pereira, E.S.; Da Silva, J.C.; Brandao, T.A.; Rocha, W.R. Phosphorane lifetime and stereo-electronic effects along the alkaline hydrolysis of phosphate esters. Phys. Chem. Chem. Phys.; 2016; 18, pp. 18255-18267. [DOI: https://dx.doi.org/10.1039/C6CP01536K]

42. Prasad, B.R.; Plotnikov, N.V.; Warshel, A. Addressing open questions about phosphate hydrolysis pathways by careful free energy mapping. J. Phys. Chem. B; 2013; 117, pp. 153-163. [DOI: https://dx.doi.org/10.1021/jp309778n]