Mechanistic insights into the phosphoryl transfer reaction in cyclin-dependent kinase 2: A QM/MM study

Cyclin-dependent kinase 2 (CDK2) is an important member of the CDK family exerting its most important function in the regulation of the cell cycle. It catalyzes the transfer of the gamma phosphate group from an ATP (adenosine triphosphate) molecule to a Serine/Threonine residue of a peptide substrate. Due to the importance of this enzyme, and protein kinases in general, a detailed understanding of the reaction mechanism is desired. Thus, in this work the phosphoryl transfer reaction catalyzed by CDK2 was revisited and studied by means of hybrid quantum mechanics/molecular mechanics (QM/MM) calculations. Our results suggest that the base-assisted mechanism is preferred over the substrate-assisted pathway when one Mg2+ is present in the active site, in agreement with a previous theoretical study. The base-assisted mechanism resulted to be dissociative, with a potential energy barrier of 14.3 kcal/mol, very close to the experimental derived value. An interesting feature of the mechanism is the proton transfer from Lys129 to the phosphoryl group at the second transition state, event that could be helping in neutralizing the charge on the phosphoryl group upon the absence of a second Mg2+ ion. Furthermore, important insights into the mechanisms in terms of bond order and charge analysis were provided. These descriptors helped to characterize the synchronicity of bond forming and breaking events, and to characterize charge transfer effects. Local interactions at the active site are key to modulate the charge distribution on the phosphoryl group and therefore alter its reactivity.


Introduction
Cyclin-dependent kinases is a family of Serine/Threonine kinases that phosphorylate peptide substrates using adenosine triphosphate (ATP) as phosphate source with a unique function in the regulation of the cell cycle [1]. As their names state, they depend on the binding of a cyclin protein in order to be fully activated [2][3][4] and also on the phosphorylation of specific residues [5][6][7][8][9]. In particular, cyclin-dependent kinase 2 (CDK2), which can bind Cyclin E or A, needs a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 obtained (42 and 24 kcal/mol, respectively) [24]. These results proposed a structural role for Asp127, however, the residue was not included in the QM region, and therefore its participation as a catalytic base could not be ruled out. More recently, Smith et al. [25] by means of QM/MM calculations, found a free energy barrier of 10.8 kcal/mol for the base-assisted mechanism and an energy barrier that could be higher than 30 kcal/mol for the substrate-assisted mechanism, since the energy barrier was not estimated once they observed a monotonically increment of the energy along the reaction pathway above 30 kcal/mol. As a conclusion, they proposed the base-assisted mechanism as the most favorable one, despite the fact that they could not estimate accurately the barrier for the substrate-assisted mechanism. It is noteworthy that the calculated energy barrier was somewhat lower than the experimentally derived value (using transition state theory), which would amount to 15.3 kcal/mol (k 3 = 35 s -1 ) [10]. At this point, they recall that this kinetic parameter was estimated using a different peptide substrate that showed a different amino acidic sequence. For instance, it had a Thr instead of a Ser residue to be phosphorylated. The authors proposed that this fact could explain the discrepancy. Among the limitations that this last study revealed is the small QM region used (only 48 atoms), which did not include most of the residues coordinating the Mg 2+ cofactor, this due to the costly sampling at a DFT level. Nonetheless, it is known that this simplification could lead to overpolarization of the QM fragments, and specially in this case the coordinating metal, affecting the energetics of the system [28]. Furthermore, the free energy barrier is dependent on the selection of the reaction coordinate, which they showed can vary considerably depending on how the combination of the reactive distances is made.
Thus, in this study the goal is to revisit the phosphoryl transfer mechanism in CDK2 with one Mg 2+ ion, allowing a more direct comparison with previous computational studies; and also to compare both mechanisms with the same QM/MM methodology, aiming at clarifying discrepancies observed in previous investigations. Besides, a detailed analysis of bond orders and atomic charges is performed in order to characterize in detail the nature of the transition states and every step of the mechanisms.

Protein preparation and molecular dynamics (MD) simulations
The initial coordinates of CDK2 were retrieved from the crystal structure with PDB (protein data bank) code 1QMZ [17]. This structure contains CDK2, a bound peptide with sequence HHASPRK, cyclin A3, an ATP molecule, and Mg 2+ as cofactor. Addition of hydrogen atoms, assignation of bond orders, and correction of bond types was done with the Protein Preparation Wizard tool included in Maestro [29][30][31]. Three residues were restored to the crystal structure: Arg297 and Leu298 on the N-terminus of CDK2, and Glu174 on the C-terminus of cyclin A3. Protonation states at physiological pH were assigned using the Propka program [32]. Titratable residues were carefully analyzed by visual inspection of hydrogen networks around them. Then, a short minimization of the structure, up to a RMSD (root-mean-square deviation) convergence criterion of 0.3 Å, was performed to eliminate any steric clashes among atoms. Parameters for all residues were derived from the OPLS-2005 force field [33].
The system was immersed in a rectangular box of water molecules using the SPC model [34] and net charge neutralization was achieved by addition of two Na + ions. A buffer region of 10 Å was imposed over all sides of the protein surface. Periodic boundary conditions were set over the x, y and z axes. The default relaxation protocol implemented in Desmond 2013 [35] was used. This consists in six steps in which the system is first energy minimized by a steepest descent algorithm applying and releasing restraints over heavy atoms. After that, four short MD simulations of 12 and 24 ps are performed retaining restraints, to finally carry out an unrestrained simulation. In order to study the viability of the substrate-assisted mechanism, a 3.4 ns production run was performed and the last snapshot was chosen for further QM/MM studies. On the other hand, to study the base-assisted mechanism, an extra simulation was carried out to bring together the H atom of the Ser residue to be phosphorylated with the oxygen atom of Asp127 to a geometry that may favor the proton transfer. Here, using a 5 ns MD simulation, the distance between these two atoms, and also the distance between the phosphorus atom of the γ-phosphate with the nucleophile oxygen atom of the Ser residue, were restrained with a force constant of 20 kcal/mol•Å 2 at equilibrium distances of 1.6 and 3.4 Å, respectively. Finally, distance restraints were removed and a 2 ns run served to extract a representative structure for QM/MM calculations. It is worth noting that to all production runs a small restraint of 0.1 kcal/mol•Å 2 was applied on the backbone atoms to keep the secondary structure of the protein stable. A 9 Å cutoff was used for the evaluation of van der Waals and electrostatic interactions, while long-range electrostatic interactions were treated with the particle mesh Ewald (PME) method [36]. Pressure was kept constant at 1 atm and temperature at 300 K using the Martyna-Tobias-Klein barostat [37] and the Nosé-Hoover chain thermostat [38], respectively. The RESPA (time-reversible reference system propagator algorithm) integrator [39] was applied with a 2 fs time step for bonded and short-range nonbonded interactions and 6 fs for long-range electrostatics.

QM/MM calculations
The last snapshot from MD simulations was used for setting QM/MM calculations. The QM region consisted in the triphosphate moiety of the ATP molecule, side chains of residues Asp145, Asn132, Lys129, Asp127, Ser5 (residue to be phosphorylated), the Mg 2+ ion, and a coordinating water molecule (Fig 2). All QM/MM cuts were done between the Cβ and Cα atoms and valences were saturated with hydrogen atoms (link-atom approach), resulting in a total of 61 atoms in the QM region. Transition state optimizations were done with a micromacro iterative procedure [40]. After that, the intrinsic reaction coordinate (IRC) path was traced down to reactant and product states, with full optimization of them. To speed-up the calculations, only the atoms within a radius of 25 Å from the ATP molecule were allowed to move, while the rest of MM atoms were kept frozen. Treatment of QM atoms was done at the B3LYP/6-31+G � [41] level of theory, similar to what has been used in previous investigations [25], while the classical OPLS-AA force field [42] was used for MM atoms, as implemented in the fdynamo library [43]. Electronic properties such as charges and Wiberg bond orders were computed for each point along the IRC path using the natural bond orbital method (NBO) [44,45]. All calculations were done with fdynamo [43] interfaced with the Gaussian 09 program [46].

Potential energy barriers
The potential energy barriers for the substrate-assisted and base-assisted mechanisms were estimated and their values are shown in Fig 3. It is possible to see how the substrate-assisted route exhibits a much higher activation energy (50.5 kcal/mol) compared to the base-assisted pathway (14.3 kcal/mol). Also, for the former mechanism, only one TS was found, making this route concerted. In the other case, two TSs were identified, which are flanked by an intermediate structure; however, the energy difference between the first transition state (TS1) and the intermediate structure (Int) is of only 1.1 kcal/mol. The energy difference between Int and TS2 is much more marked (7.2 kcal/mol), where the product state is finally reached with a reaction energy of 2.5 kcal/mol, which is very similar to the substrate-assisted mechanism (2.3 kcal/ mol). Here, and as it was mentioned, the energy difference between TS1 and Int is very small, and therefore the stepwise/concerted nature of the mechanism could change when moving from potential energy surfaces to free energy surfaces. Ongoing work in our lab is being conducted to clarify this point.
According to the calculated energy barriers, the base-assisted mechanism is more favorable than the substrate-assisted pathway when at least one Mg 2+ ion is present in the active site, what is in agreement with the last computational study carried out by Smith et al. [25]. It is worth noting that the calculated activation energy for the base-assisted mechanism (14.3 kcal/ mol) agrees very well with the experimental derived value of 15.3 kcal/mol, making this mechanism the most probable one. Regarding the substrate-assisted mechanism, our value for the activation energy (50.5 kcal/mol) could be overestimated, since other QM/MM studies in CDK2 and PKA have shown that energy barriers for this pathway are within a range of 20-40 kcal/mol [24,[47][48][49]. Here, there are also purely QM static [50] and QM/MM computational reports [48] in PKA that have shown activation energies higher than 45 kcal/mol, which would represent closer values to the ones obtained in the present study. On the other hand, an additional test was performed carrying out single point calculations on the optimized geometries using the M06-2X functional [51], but a very similar energy barrier was obtained (52 kcal/ mol). In this context, a plausible explanation for the high-energy barrier is the formation of a highly strained four-membered ring (P γ -O 1γ -H γ -O γ , see Fig 1) at the TS, effect that has also been proposed in other studies [47,50]. It is, therefore, possible that entropic effects may be especially important for modeling the substrate-assisted mechanism, and that a dynamic (QM/MM MD simulations) description of the system may be required to include the transition state flexibility, which could favor lowering the free energy barrier via entropy. This effects could, in part, explain the discrepancy with previous QM/MM calculations in CDK2 that have proposed a barrier of 24 kcal/mol for the substrate-assisted mechanism [24]. Ongoing calculations are being performed to clarify this point.

Structural analysis of reactants, TSs and products
The reactants state obtained for both reaction mechanisms was different since the conformation used was extracted from different MD simulations. In the case of the substrate-assisted mechanism, the Ser residue from the peptide substrate is establishing a hydrogen bond contact with the oxygen O 1γ of the γ-phosphate (see Figs 1 and 4A), thereby allowing the proton transfer reaction to take place to this atom. In the case of the base-assisted mechanism, an extra simulation was performed to generate a hydrogen bond between the hydroxyl group of the substrate Ser residue and residue Asp127, which showed to be stable during the MD simulation (S1 Fig). The optimized geometries of the reactants states for both reaction mechanisms are shown in Figs 4A and 5A, respectively. In the case of the substrate-assisted mechanism, the initial distance between the P γ atom and the entering oxygen O γ is 3.98 Å, slightly longer than the crystallographic distance (3.68 Å, Table 1). The TS in the substrate-assisted mechanism ( Fig 4B) shows that the phosphoryl group is midway between the donor O 3β atom and the acceptor O γ atom from the Ser residue, but pushed towards the entering oxygen atom (O 3β -P γ = 3.04 Å vs P γ -O γ = 2.10 Å, Table 1), resembling an almost planar geometry. Using Pauling's formula to assess the degree of associativity for a phosphoryl transfer reaction [52], D(n) = D(1)-0.60log(n), where D(1) is the P-O distance for a single bond (1.73 Å), and with D(n) being the average between the two P-O distances (2.57 Å), the fractional bond number (n), which gives a measure of the associativity of the mechanism, can be calculated. The value obtained for n amounts to 0.039, which means the TS is 3.9% associative or 96.1% dissociative. This result shows that though the mechanism is concerted, the TS is quite dissociative or loose, which agrees with the rather large distance between the leaving and entering oxygen atoms at the TS (4.91 Å). On the other hand, at the TS, the proton transfer reaction has also started but is not totally completed, showing the H atom between both donor and acceptor oxygen atoms ( Fig 4B). Fig 4C shows the product state, where it is possible to observe how the phosphoryl transfer reaction has been completed, with the characteristic that the transferred proton is forming a hydrogen bond with one of the oxygen atoms of the β-phosphate (O 3β ) of the newly formed ADP molecule. The Lys129 and Asp127 residues help in maintaining the in-line configuration necessary for the phosphoryl transfer reaction to take place and, especially, Lys129 stabilizes the transferred phosphoryl group at the TS through hydrogen bonding interactions. It is also possible to observe that the octahedral coordination around the Mg 2+ ion is preserved during and after the reaction (Table 1), in agreement with previous computational studies [24,25].
In the case of the base-assisted mechanism, Fig 5B shows how in the first TS (TS1) the phosphoryl group has a planar geometry with almost equal distances between P γ and the leaving and entering oxygen atoms (2.57 and 2.47 Å, respectively). The distance between those two oxygen atoms is 5.03 Å, revealing the high dissociative character of the mechanism, and with a fractional bond number of 0.048, i.e., 95.2% dissociative. Very close in energy is the intermediate structure (Int, Fig 5C), which also resembles a metaphosphate geometry but with a shorter distance to the entering oxygen atom (2.12 Å, Table 1). At this point, the distance of the oxygen atom O 2γ coordinating the Mg 2+ ion has increased from 2.06 Å to 2.39 Å. The reaction proceeds and a second transition state (TS2) is reached, where the P γ -O γ bond is almost formed (1.83 Å). At this stage, the proton transfer to Asp127 has been initiated and the proton is midway between the donor (O γ ) and acceptor (O δ1 ) oxygen atoms. An interesting feature of this Phosphoryl transfer reaction in cyclin-dependent kinase 2: A QM/MM study structure is the spontaneous transfer of a proton from Lys129 to the O 2γ atom of the ATP phosphoryl group. This proton transfer reaction is almost complete at this transition state, with the proton lying between the donor nitrogen and acceptor oxygen atoms ( Fig 5D). Here, it is very probable that neutralization of Asp127 plays a role in lowering the pK a of Lys129, favoring in this way the proton transfer. Finally, the system reaches the product state (Fig 5E), which features the phosphorylated Ser residue with a distance between P γ to the O 3β oxygen of 3.66 Å. Asp127 is protonated forming a hydrogen bond with the oxygen atom of the Ser residue. The phosphorylated Ser residue is now protonated at the O 2γ oxygen and forming a hydrogen bond with the nitrogen atom of Lys129. The O 2γ oxygen has partially lost its coordination with the Mg 2+ ion, since the distance to it has increased to 3.18 Å. All the rest of the  (Table 1). Finally, it is expected that the usual protonation states of Asp127 and Lys129 may be easily recovered when the products leave the active site.
These results show that while the substrate-assisted mechanism exhibits a concerted character, the base-assisted mechanism is dissociative. Asp127 serves as a base accepting the proton from the substrate Ser residue late in the progress of the reaction, what also agrees with the last computational study in CDK2 [25] and with QM/MM calculations in PKA [47]. It is interesting that in both mechanisms the TSs exhibit a dissociative character, though the substrateassisted mechanism has been related to more associative mechanisms [24,47]. The possible cause for this could be the absence of the second Mg 2+ in our simulations, which would reduce the repulsion between the reactive fragments. On the other hand, one of the new features found in the base-assisted mechanism is the proton transfer from Lys129 to the transferred phosphoryl group. This result would suggest that Lys129 not only exerts its structural function bridging closer the γ-phosphate and the Ser residue, but could also help to stabilize the negative charge created on the phosphoryl group at the TS by a proton transfer reaction. A similar charge stabilizing event was also observed by QM/MM calculations in Glucokinase [53], where Lys169 was identified as an acid catalyst protonating one of the oxygen atoms of the γphosphate in the ATP molecule. The great structural similarity among protein kinases allows for some comparison of the mechanisms since both enzymes share the common residues in the active site that are key for enzyme catalysis. In the case of CDK2, as mentioned previously, it is now postulated that two Mg 2+ ions are required for optimal catalysis [21,22]. Despite that this structural effect is out of the scope of the present investigation, it is plausible to think that in the absence of a second Mg 2+ ion, which would neutralize the negative charge in the active site, a proton transfer reaction from Lys129 may help to mimic its charge stabilizing function.
In order to complement the analysis of the base-assisted mechanism, we have performed QM calculations to predict the pK a values of Asp127, Lys129, and the phosphorylated serine group (pSer) at the O 2γ atom using the product complex structure. The complete QM cluster model consisted in the Mg 2+ ion with its coordinating water, the ADP molecule, the Phosphoryl transfer reaction in cyclin-dependent kinase 2: A QM/MM study phosphorylated serine residue (pSer), and residues Asp127, Lys129, Asn132, Asp145, and Thr165. We performed different types of calculations to test how the inclusion of different residues in the QM model affected the calculated pK a 's, and all the results have been summarized in S1 Table. For instance, it can be observed how the experimental pK a values of Lys, Asp and phosphorylated serine (pSer) residues in water are well reproduced with the adopted protocol (10.2, 3.5 and 5.8, respectively.). Since we were interested in computing pK a values at specific geometries, different constraints were used that are described in the supplementary information (S1 Protocol). The pK a for the protonated oxygen atom of Asp127 was calculated to be 6.8 in the complete model, what reveals that the active site of CDK2, and in particular its interaction with the pSer group, produces a shift towards higher values; removing the pSer group from the QM model restores a pK a of 3.9. This shift is expected in terms of electrostatic interactions since Asp127 is now closer to the transferred γ-phosphate group, and Lys129 is in its neutral state. On the other hand, we evaluated the pK a value of the pSer group (at the oxygen atom O 2γ ) with the surrounding active site residues except for Lys129. This calculation results in a pK a value of 14.5, showing a significant shift with respect to the pK a value in water (5.8).
Here, the presence of negatively charged groups like ADP and Asp145 may help to generate such a high pK a value and would explain why the phosphate oxygen is a suitable proton acceptor in the base-assisted mechanism. Since Lys129 is in a deprotonated state in the product complex, the inclusion of this residue further increases the phosphate oxygen's pK a to 20.2, since a hydrogen bond, with the nitrogen atom of Lys129, stabilizes the proton in O 2γ . Thus, what these calculations suggest is that the γ-phosphate oxygen is a strong base ready to accept a proton from the surroundings. The pK a in Lys129 was also evaluated in a model where the pSer fragment was not included, which resulted in a very low pK a (4.3). This shift could be explained by the direct microenvironment around Lys129 (Asn132, Asp127, and Thr165), which is quite neutral in the product state, also considering the positively charged Mg 2+ ion, and with the presence of Asp145 and ADP in a more distant position. The inclusion of the pSer fragment, which is in a protonated state, further diminishes Lys129's pK a at the product conformation to a value of 1.1. It is important to emphasize that this very low value is obtained because the residues' geometries are restricted, and only relaxation of the NH 3 group of Lys129 is allowed. Therefore, Lys129 in that particular conformation is highly stabilized in a deprotonated state, forming a hydrogen bond with the protonated oxygen O 2γ . A pK a = 10.9 was obtained under two conditions: first, the entire model system relaxation and second, the optimization of phosphate group at pSer residue, with Lys129 in its protonated and deprotonated states. This new value comes due to the rotation of the phosphate group in the protonated state of Lys129, where the protonated phosphate oxygen interacts with the β-phosphate oxygens of ADP and protonated Lys129 interacts with one deprotonated oxygen at γ-phosphate. This new interaction stabilizes the protonated state of Lys129 producing the calculated pK a . This also shows that just upon rotation of the phosphate group, Lys129 could quickly restore its original protonation state since in the presence of a negatively charged phosphate oxygen there is a marked shift in its pK a . A possibility could be a proton transfer from Asp127 or a water molecule in the active site. We want to stress that these results should not be considered as exhaustive or taken in a quantitative sense. The product complex was chosen since is structurally similar to the second transition TS2 where the proton transfer from Lys129 takes place. However, the pK a of the different residues in the active site may change when going from reactants to products. However, our results would be suggesting that the γ-phosphate oxygen O 2γ exhibits a high pK a that can favor the proton transfer from Lys129. The approach of the γ-phosphate to the serine residue would also increase the pK a of Asp127, making it a suitable proton acceptor. In this context, theoretical investigations have also proposed that high pK a values (from 11 to 15) [54][55][56] in non-bridging oxygens of phosphorane intermediates in ribozyme catalysis may promote proton transfers, resembling a similar scenario to the one postulated in the present investigation. Finally, we want also to point out that the plausibility of the proposed mechanism should be confirmed by QM/MM free energy calculations since our methodology largely depends on the initial conformation that was extracted from classical MD simulations, which positioned Lys129 very close to the γ-phosphate of ATP, which also could have favored the proton transfer from this residue.

Synchronicity studied by bond order analysis
As previously described, single point calculations using the NBO method [57] on each point along the reaction paths were performed in order to calculate Wiberg bond orders. Bond orders are a direct measurement of the formation or breakage's degree of a bond and therefore are a good property to describe the synchronicity of bond breaking and bond forming events. In order to get clearer insights into these events, it has been postulated the use of bond order derivatives, as a useful qualitative property to describe synchronicity and electronic changes in chemical reactions [58,59]. These will show positive values for bond formation or strengthening, while negative values are expected for bond breaking or weakening. Fig 6A and 6B show the evolution of the bond orders for both mechanisms along the reaction coordinate, while Fig  6C and 6D show bond order derivatives. It is readily seen from Fig 6A that bond orders follow the expected behavior for the substrate-assisted reaction: the bond O 3β -P γ begins to break gradually while the rest of the bonds are maintained almost unaltered until this bond is completely broken; then, the bond P γ -O γ begins to form. This supports the fact that P γ at the TS structure shows a longer distance with respect to the leaving oxygen (O 3β ) than with the entering oxygen (O γ ). At the TS, the curves representing the bonds O γ -H γ , O 1γ -H γ and P γ -O γ cross with a value for the bond order of approximately 0.6 (S2 Table), showing that at the TS the formation and cleavage of these bonds has advanced halfway. This also shows that there is a marked synchronicity in the formation and dissociation of these bonds. However, this is better illustrated in Fig 6C, where bond order derivatives show an initial negative sign for the O 3β -P γ bond, depicting its gradual breaking, and how the peaks for the other three bonds are located exactly at the TS, with the positive peaks representing the formation of the O 1γ -H γ and P γ -O γ bonds, and a negative peak representing the breaking of the O γ -H γ bond. Thus, these results show that though the reaction is concerted, there is a high degree of asynchronicity with respect to the breaking of the O 3β -P γ bond, while the rest of the events occur in a synchronous fashion. Furthermore, it is seen how the peak representing the formation of the P γ -O γ bond is quite broad in direction to the product state, showing how the strengthening of this bond is gradual. Finally, the last negative peak in the O 1γ -H γ bond represents some elongation that this bond undergoes when start forming a hydrogen bond with the O 3β oxygen of ADP.
In the case of the base-assisted mechanism, bond breaking and forming events occur in a similar way, beginning with the breaking of the O 3β -P γ bond (Fig 6B) until TS1 is reached. At this stage, the bond P γ -O γ has just started to form and the reaction proceeds with the strengthening of this bond until the intermediate structure is formed. At this point, this bond is almost half-formed, i.e., bond order of 0.46 (S3 Table), while the bond O 3β -P γ is completely broken. Following the reaction coordinate, it is possible to observe small variations in the bond orders for the rest of the bonds due to structural rearrangements prior the completion of the phosphoryl and proton transfer reactions. As was discussed previously, at TS2 the phosphoryl transfer is almost complete, with a bond order value of 0.83 for the P γ -O γ bond and with intermediate values for the rest of the bonds that participate in both proton transfers (S3 Table). Fig 6D shows clearly the degree of synchronicity of bond forming and breaking events. As expected, the first negative peak represents the breaking of the O 3β -P γ bond, followed by positive values representing some degree of formation of the P γ -O γ bond. In this case, bond formation with O γ occurs markedly before reaching TS2, where the proton transfer reactions take place. It is also observed that the proton transfer from Lys129 to the phosphoryl group occurs in a slightly asynchronous way with respect to the proton transfer to Asp127 and the strengthening of the P γ -O γ bond, i.e., though the peaks are very close, they do not have their maxima at the same position on the reaction coordinate. The positive and negative peaks representing the proton transfer from Lys129 to the phosphoryl group are located slightly before the peaks representing the proton transfer to Asp127 and the completion of the P γ -O γ bond. This is also  Table).

Charge analysis
In order to gain a deeper insight into the differences between the two explored reaction paths and to evaluate possible charge transfer effects, NPA (natural population analysis) charges [60] were calculated and are plotted in Fig 7. Fig 7A and 7B show the evolution of the charge on the P γ atom along the reaction coordinate. It is observed that the trend of the atomic charge on this atom follows the same behavior in both mechanisms: the charge decreases (more negative values) as the mechanism approaches the TS region and increases its value until the product state. In the case of the substrate-assisted mechanism, the initial value is 2.55 (charges are presented in a.u. in this section) which decreases to a value of 2.49 before reaching the TS. At the TS, it takes a value of 2.52, from where increases until reaching a value of 2.61 at the product state (S4 Table). Thus, though there are changes on the charge of this atom, these fluctuations are not large. This event is accompanied by a shallow increase in the charge of the Mg 2+ ion (Fig 7E), showing that this ion loses some electron density as the coordination with O 2γ is weakened. It is also seen that the charge on the non-bridging oxygen atoms increases (more positive) slightly and then decreases before reaching the TS, especially for O 1γ and O 3γ ( Fig  7C). This behavior can be explained since the breaking of the O 3β -P γ bond withdraws electron density from the P γ atom and charge transfer from the non-bridging oxygen atoms would take place to compensate for the created electron deficiency in P γ . On the other hand, it is also observed that the non-bridging oxygens O 2γ and O 3γ possess a similar atomic charge, but with O 2γ bearing a slightly more negative charge, since it is expected to be more polarized by the coordination with the Mg 2+ ion. On the other hand, the oxygen atom O 1γ has a more positive charge, most probably since it is only stabilized by one hydrogen bond with the hydroxyl group of the Ser residue ( Fig 8A). Otherwise, O 3γ is initially stabilized by three hydrogen bonds established with three water molecules, providing a larger stabilization for the charge developed on this atom (Fig 8A). However, the hydrogen bond that is initially formed with the water molecule that coordinates Mg 2+ is lost before reaching the TS. It is worth to note that variations in the atomic charges are more pronounced on the oxygen atoms compared to the changes on the P γ atom and the Mg 2+ ion, and therefore the changes on the phosphoryl group (PO 3 -) charge will be determined by the charge fluctuations on the non-bridging oxygens. Before reaching the TS, the charges on O 2γ and O 3γ are the same, since at this stage the phosphoryl group resembles a metaphosphate-like structure, this is, the bond order O 3β -P γ is close to zero and the new P γ -O γ bond has not yet been formed.
With these observations at hand, it becomes clear that the effect that the protein environment exerts on the charge distribution of the phosphoryl group is very important and that charge stabilization depends on specific interactions at the active site, especially on hydrogen bonding interactions. At the TS, the charge on the phosphoryl group increases (more positive) sharply due to the proton transfer to the O 1γ atom. After the TS, the charge on the phosphoryl group remains stable and decreases before reaching the product state, accounting for the formation of the P γ -O γ bond. The charges on O 2γ and O 3γ also become more negative due to the formation of the new bond, but now the charge on O 2γ becomes slightly more positive than O 3γ , most probably since the distance with Mg 2+ has increased and hydrogen bond stabilization on O 3γ is enhanced. It is also observed that the charge on Mg 2+ increases monotonically (Fig 7E), accounting for the continuing weakening of its coordination with O 2γ , until the formation of the product state; except for a small fluctuation at the TS, due to the change in the charge distribution of the phosphoryl group upon proton transfer. In the case of the base-assisted mechanism, the P γ atom also decreases its charge (less positive) as approaches the TS (Fig 7B). The initial value is 2.61 and then decreases to 2.57 at TS1, and it diminishes slightly more at the intermediate structure with a value of 2.56 (S4 Table). At TS2 the value has increased to 2.59, with a final value at the product state of 2.60. Interestingly, compared to the substrate-assisted mechanism, the difference in the charge value between the reactant state and the lowest value along the reaction coordinate is the same (0.06); however, the magnitude of the charge on this atom is more positive in the second path, showing that the P γ atom in the base-assisted mechanism is more reactive to the attack of the entering oxygen atom, since its electrophilicity is enhanced. This effect could also help to explain the lower energy barrier for the base-assisted mechanism. As with the previous analysis, the charge on the Mg 2+ ion smoothly increases along the reaction coordinate (Fig 7F), but in this case the difference between the charges, at the reactant and product state, is larger (0.07) when compared to the substrate-assisted mechanism (0.04). On the other hand, it is also seen that the charge on the non-bridging oxygens is slightly more negative compared to the substrate-assisted mechanism, what also makes the charge on the phosphoryl group more negative (Fig 7D). Specifically, the charge on O 3γ is notoriously more negative when compared to O 1γ and O 2γ . This is explained since O 3γ is stabilized by three hydrogen bonds from three water molecules (Fig 8B), what brings increased charge stabilization for that atom. In the case of O 1γ , upon the absence of a hydrogen bond with the Ser residue, two hydrogen bonds with two water molecules are formed instead ( Fig 8B); this results in a charge slightly more positive than for the O 2γ atom, which coordinates to the Mg 2+ ion. Hence, the charge on the phosphoryl group follows the trend of the non-bridging oxygen atoms; this is, the charge increases (more positive) until reaching TS1, where it takes a value of -0.98. As discussed previously, this effect is due to the dissociation of the O 3β -P γ bond, event that withdraws electron density from the phosphoryl group. After that, the charge begins to gradually decrease until reaching the intermediate structure with a charge of -1.03. This is due to the fact that some degree of formation of the P γ -O γ bond is already observed at this stage (see Fig  6D and its discussion), what generates a gradual recovery of the negative charge on the phosphoryl group. From this point, the charge value is practically maintained until just before reaching TS2, what agrees with the little formation of the P γ -O γ bond between Int and TS2. Then, the charge increases to -0.94 due to the protonation of O 2γ by Lys129, and then diminishes slightly to a value of -0.97 at the product state, due to the completion of the P γ -O γ bond. https://doi.org/10.1371/journal.pone.0215793.g008 Phosphoryl transfer reaction in cyclin-dependent kinase 2: A QM/MM study It is interesting to note that the charge on the phosphoryl group at TS1 and Int is close to the charge expected in a dissociative mechanism: -1 for a metaphosphate species [20]. The more negative charge on the phosphoryl group in the base-assisted mechanism correlates well with a larger change in the charge on Mg 2+ (Fig 7F), which donates electron density to the phosphoryl group upon dissociation of the coordinating bond with O 2γ .
In this way, our results highlight the importance of local interactions in the modulation of charge on the phosphoryl group. Furthermore, it was shown how in a reaction with high dissociative character the charge on the phosphoryl group becomes more positive when it reaches a metaphosphate species, an issue that was not clearly established in the previous literature [20]. On the contrary, for an associative mechanism, a negative charge buildup at the TS is expected, a fact that has been previously characterized computationally [61]. Our results support experimental observations using 19 F NMR on MgF 3 transition state analogues, which have been able to describe how the direct microenvironment in near transition-state conformations alter the charge distribution of the TS mimic, i.e., by local electrostatic and hydrogen bonding interactions [62]. This is exactly what our results show for the phosphoryl group in both mechanisms: the charge distribution is finely modulated by local interactions with the Mg 2+ ion, Lys129, the substrate Ser residue, and a network of water molecules. This aspect is expected to be crucial for the catalytic mechanism of protein kinases since charge stabilization at the TS is recognized to be a critical point for enzyme catalysis. In particular, the crystal structure of CDK2 used in this study contains the glycine-rich (Gly-rich) loop in the open state; this structural moiety functions as a lid closing the active site. It has been observed that the presence of two Mg 2+ ions in the active site favors the closed state [21], where backbone amides of the Gly residues make hydrogen bonding interactions with the β-phosphate and presumably with the phosphoryl group at the TS. Thus, it is expected that in the presence of two Mg 2+ ions the charge distribution on the phosphoryl group may be stabilized by specific hydrogen bonding and electrostatic interactions that favor the catalysis. However, how and how much these charge stabilizing interactions affect the phosphoryl transfer reaction is a question that remains open. These points are expected to be covered in future work.

Conclusions
CDK2 is a very important kinase that is considered a study model for the cyclin-dependent kinase family. As other kinases, this enzyme catalyzes the phosphoryl transfer reaction from an ATP molecule to a peptide substrate containing a Ser o Thr residue to be phosphorylated.
Though this system has been subject of experimental and computational studies, the two proposed reaction mechanisms, substrate-assisted and base-assisted, had not been studied until now within the same computational approach, which is a requisite in order to have a proper comparison between them. In this context, our results suggest that the base-assisted mechanism is the preferred path, at least when one Mg 2+ ion is used as cofactor, with an estimated potential energy barrier of 14.3 kcal/mol, which is in good agreement with the experimentally derived value. This mechanism is stepwise, in contrast to the substrate-assisted mechanism that is concerted. Both mechanisms show TSs that have a high dissociative character. Interestingly, a new feature in the base-assisted mechanism has been observed: a spontaneous proton transfer from Lys129 to one of the oxygen atoms of the transferred phosphoryl group. This event takes place late in the reaction progress, at the second TS, with the proposed effect of neutralizing the negative charge on the phosphoryl group. We recognize that a major limitation in the present work is the absence of a second Mg 2+ ion in the active site, which is expected to modulate the energy barrier and possibly the mechanism. These aspects are expected to be covered in future work, and therefore, the prevalence of the base-assisted mechanism remains to be confirmed. Despite this, the calculated energy barrier seems to be reasonable, and agrees with the experimental data, what would be suggesting that the reaction could still be carried out with only one Mg 2+ ion in the active site. Detailed mechanistic information has been given in terms of bond orders and natural population analysis. In both mechanisms, the breaking of the O 3β -P γ bond precedes the formation of the new bond with the entering oxygen from the substrate Ser residue. The rest of the bond breaking and bond forming events, related with the proton transfer to the phosphate group or to residue Asp127, occur in a synchronous fashion, although some degree of asynchronicity is observed for the base-assisted mechanism. Charge analysis showed that, due to the dissociative character of the mechanisms, the charge on the phosphoryl group becomes more positive when a metaphosphate species is formed. These results help to clarify how is the charge evolution on key atoms during the enzyme catalysis, such as the Mg 2+ ion, the P γ atom and the nonbridging oxygens of the phosphoryl group. Charge transfer effects appear to have an important role, specifically on the Mg 2+ ion, which delivers electron density as the reaction moves to the product state. Also, it was observed that atomic charges on the non-bridging oxygens are modulated by local interactions from the protein microenvironment, mainly electrostatic and hydrogen bonding interactions with, primarily, the Mg 2+ cofactor, active site residues such as Lys129, and water molecules. Future work will be concentrated in elucidating the effect of a second Mg 2+ ion in the mechanism and its effect on the charge distribution of the phosphoryl group. In summary, these results help to advance towards a more comprehensive understanding of phosphoryl transfer reactions in protein kinases, information that can be highly relevant for the development of new inhibitors.  Table. Calculated pK a values for Lys129, Asp127 and pSer residues in the product state in the presence and absence of different active site residues. a reference [63]. b reference [64]. (DOCX) S2 Table. Calculated Wiberg bond orders, for the bonds that are broken and formed, for the substrate-assisted mechanism at the three stationary points. (DOCX) S3 Table. Calculated Wiberg bond orders, for the bonds that are broken and formed, for the base-assisted mechanism at the five stationary points. (DOCX) S4 Table. NPA charges (e) for the most relevant atoms and the PO 3 group for both mechanisms at the different stationary points. (DOCX) S1 Protocol. Protocol adopted for the calculation of pK a values of different residues of the active site in a QM cluster model. (DOCX)