Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Open Access

Peer-reviewed

Research Article

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

Abstract

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.

Citation: Recabarren R, Osorio EH, Caballero J, Tuñón I, Alzate-Morales JH (2019) Mechanistic insights into the phosphoryl transfer reaction in cyclin-dependent kinase 2: A QM/MM study. PLoS ONE 14(9): e0215793. https://doi.org/10.1371/journal.pone.0215793

Editor: Matthias Stein, Max Planck Institute for Dynamics of Complex Technical Systems, GERMANY

Received: April 5, 2019; Accepted: August 19, 2019; Published: September 4, 2019

Copyright: © 2019 Recabarren et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the paper and its Supporting Information files.

Funding: J.A.M. acknowledges the financial support from FONDECYT projects 1140618 and 1181253. R.R. acknowledges support from the doctoral fellowship CONICYTPFCHA/Folio 21130949. I.T. thanks the financial support from Generalitat Valenciana through project AICO/2018/238. The authors also acknowledge to “Centro de Bioinformática, Simulación y Modelado (CBSM)” at Universidad de Talca for providing computational resources to carry out the calculations reported in this study.

Competing interests: The authors have declared that no competing interests exist.

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–4] and also on the phosphorylation of specific residues [5–9]. In particular, cyclin-dependent kinase 2 (CDK2), which can bind Cyclin E or A, needs to be phosphorylated at Thr160 [10]. CDK2 helps in the progression from G1 to S phase during the cell cycle and its malfunctioning, e.g. by mutations, has been related to different human cancers [11]. For this reason, many CDK inhibitors have been proposed, which in the majority of the cases bind in the ATP binding site. Also, and most recently, new inhibitors have been designed to tackle protein-protein interactions and allosteric sites [12]; however, the development of potent and selective inhibitors has been very challenging and in many cases disappointing, without getting through clinical trials [13]. In this context, the precise knowledge of the phosphoryl transfer mechanism in CDKs, and kinases in general, could lead to new strategies for the development of more potent and selective drugs. Moreover, CDK2 is a very interesting system to be used as a model of study within the CDK family due to the availability of many crystallographic structures and kinetic data [10,14–18].

In general, CDKs and kinases provide specific amino acids and Mg2+ cofactors that help to position ATP and the substrate to be phosphorylated in the correct orientation for catalysis. Also, the active site provides compensatory charges that stabilize the negative charge at the transition state (TS), lowering in this way the energy barrier of the reaction [19,20]. However, many aspects of the phosphoryl transfer mechanism are still unclear, such as the role of Mg2+ ions in the mechanism and which is the most favorable catalytic pathway [20]. In the former aspect, CDK2 was always thought to function efficiently with only one Mg2+ ion within the active site, different to what is known for other kinases, such as PKA (protein kinase A), that uses two Mg2+ ions as cofactors. However, recent crystallographic and kinetic data have shown that CDK2 also works more efficiently with two Mg2+ ions. This second Mg2+ ion prevents ADP release and therefore it needs to be removed from the active site after the phosphoryl transfer step is completed [21,22].

On the other hand, when it comes to the phosphoryl transfer mechanism itself, two main scenarios are possible: one in which the proton from the hydroxyl group of the Ser/Thr side chain is transferred to a nearby aspartate residue (Asp127 in CDK2), route that is called base-assisted mechanism; or a second one where the proton can be transferred to one of the oxygen atoms of the γ-phosphate of the ATP molecule, pathway that is called substrate-assisted mechanism (Fig 1) [23–25]. These two mechanisms are associated to another classification of phosphoryl transfer mechanisms, named as dissociative, concerted and associative mechanisms [20]. In a dissociative mechanism, the breaking of the P-O bond with the leaving oxygen atom precedes bond formation and the reaction goes through a metaphosphate intermediate linked by two TSs. In this case, the TS is dissociative in nature or “loose”, characterized by long P-O distances with respect to the leaving and entering oxygens. In an associative mechanism, the distances of the breaking and forming P-O bonds are much shorter and the reaction goes through a pentavalent phosphorane intermediate; in this way the TS is associative or “tight”. Also, the mechanism can be concerted, which implies a single TS, where the TS may have a more dissociative or associative nature depending on the bond orders of the bonds to be broken and formed. It is worth noting that a third option for the deprotonation of the Ser/Thr hydroxyl group could be a proton shuttle from water-mediated interactions as has been proposed in phosphatase reactions [26]; this option is not covered in the present investigation and therefore cannot be ruled out as a possible reaction mechanism [24].

thumbnail
Download:
Fig 1. Proposed pathways for the substrate-assisted and base-assisted mechanisms.

Dotted lines represent bonds that are broken or formed during the reaction.

https://doi.org/10.1371/journal.pone.0215793.g001

Previous computational studies in CDK2 have not completely defined which is the most favorable route [24,27]; however, insights that the base-assisted mechanism is preferred have been given more recently [25]. The first computational study used cluster calculations at a DFT (density functional theory) level and proposed the substrate-assisted mechanism as the operating one [27]. A subsequent study carried out by the same authors, but now considering a QM/MM (quantum mechanics/molecular mechanics) partition of the system, also reaffirmed that fact, with an energy barrier for the reaction much lower than the previous one 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 (k3 = 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 Mg2+ 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 Mg2+ 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.

Methods

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 Mg2+ 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–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 Mg2+ 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 micro-macro 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].

thumbnail
Download:
Fig 2. Atoms included in the QM region for QM/MM calculations.

The valences of the atoms lying at the QM/MM boundary were completed with hydrogen atoms.

https://doi.org/10.1371/journal.pone.0215793.g002

Results and discussion

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.

thumbnail
Download:
Fig 3. Potential energy profiles for both, substrate-assisted (black) and base-assisted (red) mechanisms.

The numbers show the relative energy for each stationary point with respect to reactants.

https://doi.org/10.1371/journal.pone.0215793.g003

According to the calculated energy barriers, the base-assisted mechanism is more favorable than the substrate-assisted pathway when at least one Mg2+ 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–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γ-O1γ-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 O1γ 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 O3β atom and the acceptor Oγ atom from the Ser residue, but pushed towards the entering oxygen atom (O3β-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 (O3β) 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 Mg2+ ion is preserved during and after the reaction (Table 1), in agreement with previous computational studies [24,25].

thumbnail
Download:
Table 1. Measured distances between key atoms that participate during catalysis and those that complete the octahedral coordination around the Mg2+ ion.

https://doi.org/10.1371/journal.pone.0215793.t001

thumbnail
Download:
Fig 4. Molecular structures of the stationary points in the substrate-assisted mechanism.

(A) Reactant state. (B) Transition state. (C) Product state.

https://doi.org/10.1371/journal.pone.0215793.g004

thumbnail
Download:
Fig 5. Molecular structures of the stationary points in the base-assisted mechanism.

(A) Reactant state. (B) Transition state TS1. (C) Intermediate state Int. (D) Transition state TS2. (E) Product state.

https://doi.org/10.1371/journal.pone.0215793.g005

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 O2γ coordinating the Mg2+ 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 structure is the spontaneous transfer of a proton from Lys129 to the O2γ 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 pKa 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 O3β 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 O2γ oxygen and forming a hydrogen bond with the nitrogen atom of Lys129. The O2γ oxygen has partially lost its coordination with the Mg2+ ion, since the distance to it has increased to 3.18 Å. All the rest of the O-Mg2+ distances are in general conserved in values that do not go beyond 2.2 Å (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 substrate-assisted mechanism has been related to more associative mechanisms [24,47]. The possible cause for this could be the absence of the second Mg2+ 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 Mg2+ 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 Mg2+ 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 pKa values of Asp127, Lys129, and the phosphorylated serine group (pSer) at the O2γ atom using the product complex structure. The complete QM cluster model consisted in the Mg2+ ion with its coordinating water, the ADP molecule, the 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 pKa’s, and all the results have been summarized in S1 Table. For instance, it can be observed how the experimental pKa 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 pKa values at specific geometries, different constraints were used that are described in the supplementary information (S1 Protocol). The pKa 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 pKa 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 pKa value of the pSer group (at the oxygen atom O2γ) with the surrounding active site residues except for Lys129. This calculation results in a pKa value of 14.5, showing a significant shift with respect to the pKa value in water (5.8). Here, the presence of negatively charged groups like ADP and Asp145 may help to generate such a high pKa 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 pKa to 20.2, since a hydrogen bond, with the nitrogen atom of Lys129, stabilizes the proton in O2γ. Thus, what these calculations suggest is that the γ-phosphate oxygen is a strong base ready to accept a proton from the surroundings. The pKa in Lys129 was also evaluated in a model where the pSer fragment was not included, which resulted in a very low pKa (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 Mg2+ 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 pKa 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 NH3 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 O2γ. A pKa = 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 pKa. 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 pKa. 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 pKa 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 O2γ exhibits a high pKa that can favor the proton transfer from Lys129. The approach of the γ-phosphate to the serine residue would also increase the pKa of Asp127, making it a suitable proton acceptor. In this context, theoretical investigations have also proposed that high pKa values (from 11 to 15) [54–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 O3β-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 (O3β) than with the entering oxygen (Oγ). At the TS, the curves representing the bonds Oγ-Hγ, O1γ-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 O3β-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 O1γ-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 O3β-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 O1γ-Hγ bond represents some elongation that this bond undergoes when start forming a hydrogen bond with the O3β oxygen of ADP.

thumbnail
Download:
Fig 6. Wiberg bond orders and bond order derivatives of the bonds that are broken or formed during the reaction for both mechanisms.

(A) Bond orders for the substrate-assisted mechanism. (B) Bond orders for the base-assisted mechanism. (C) Bond order derivatives for the substrate-assisted mechanism. (D) Bond order derivatives for the base-assisted mechanism. Dotted lines show the position of the stationary points on the reaction coordinate.

https://doi.org/10.1371/journal.pone.0215793.g006

In the case of the base-assisted mechanism, bond breaking and forming events occur in a similar way, beginning with the breaking of the O3β-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 O3β-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 O3β-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 revealed by the bond order values of the bonds that are formed at TS2 (O2γ-Hζ1 = 0.87 vs Oδ1-Hγ = 0.62, S3 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 Mg2+ ion (Fig 7E), showing that this ion loses some electron density as the coordination with O2γ 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 O1γ and O3γ (Fig 7C). This behavior can be explained since the breaking of the O3β-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 O2γ and O3γ possess a similar atomic charge, but with O2γ bearing a slightly more negative charge, since it is expected to be more polarized by the coordination with the Mg2+ ion. On the other hand, the oxygen atom O1γ 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, O3γ 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 Mg2+ 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 Mg2+ ion, and therefore the changes on the phosphoryl group (PO3-) charge will be determined by the charge fluctuations on the non-bridging oxygens. Before reaching the TS, the charges on O2γ and O3γ are the same, since at this stage the phosphoryl group resembles a metaphosphate-like structure, this is, the bond order O3β-Pγ is close to zero and the new Pγ-Oγ bond has not yet been formed.

thumbnail
Download:
Fig 7. Evolution of NPA charges for the most relevant atoms in the reaction.

(A, C and E) Charges for Pγ, PO3 and Mg2+ in the substrate-assisted mechanism, respectively. (B, D and F) Charges for Pγ, PO3 and Mg2+ in the base-assisted mechanism, respectively. Dotted lines show the position of the stationary points on the reaction coordinate.

https://doi.org/10.1371/journal.pone.0215793.g007

thumbnail
Download:
Fig 8. Reactants structures showing hydrogen bonding interactions (dashed blue lines) that stabilize the γ-phosphate of the ATP molecule.

(A) Substrate-assisted mechanism. (B) Base-assisted mechanism.

https://doi.org/10.1371/journal.pone.0215793.g008

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 O1γ 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 O2γ and O3γ also become more negative due to the formation of the new bond, but now the charge on O2γ becomes slightly more positive than O3γ, most probably since the distance with Mg2+ has increased and hydrogen bond stabilization on O3γ is enhanced. It is also observed that the charge on Mg2+ increases monotonically (Fig 7E), accounting for the continuing weakening of its coordination with O2γ, 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 Mg2+ 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 O3γ is notoriously more negative when compared to O1γ and O2γ. This is explained since O3γ is stabilized by three hydrogen bonds from three water molecules (Fig 8B), what brings increased charge stabilization for that atom. In the case of O1γ, 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 O2γ atom, which coordinates to the Mg2+ 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 O3β-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 O2γ 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. 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 Mg2+ (Fig 7F), which donates electron density to the phosphoryl group upon dissociation of the coordinating bond with O2γ.

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 19F NMR on MgF3- 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 Mg2+ 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 Mg2+ 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 Mg2+ 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 Mg2+ 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 Mg2+ 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 Mg2+ 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 O3β-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 Mg2+ ion, the Pγ atom and the non-bridging oxygens of the phosphoryl group. Charge transfer effects appear to have an important role, specifically on the Mg2+ 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 Mg2+ cofactor, active site residues such as Lys129, and water molecules. Future work will be concentrated in elucidating the effect of a second Mg2+ 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.

Supporting information

S1 Fig. Distance between the Hγ atom from the substrate Ser residue and the Oδ1 atom from Asp127 showing a stable hydrogen bonding interaction during an unrestrained 2 ns MD simulation.

https://doi.org/10.1371/journal.pone.0215793.s001

(PNG)

S1 Table. Calculated pKa 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].

https://doi.org/10.1371/journal.pone.0215793.s002

(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.

https://doi.org/10.1371/journal.pone.0215793.s003

(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.

https://doi.org/10.1371/journal.pone.0215793.s004

(DOCX)

S4 Table. NPA charges (e) for the most relevant atoms and the PO3 group for both mechanisms at the different stationary points.

https://doi.org/10.1371/journal.pone.0215793.s005

(DOCX)

S1 Protocol. Protocol adopted for the calculation of pKa values of different residues of the active site in a QM cluster model.

https://doi.org/10.1371/journal.pone.0215793.s006

(DOCX)

Acknowledgments

The authors also acknowledge to “Centro de Bioinformática, Simulación y Modelado (CBSM)” at Universidad de Talca for providing computational resources to carry out the calculations reported in this study.

References

  1. 1. Morgan DO. Principles of CDK regulation. Nature 1995;374:131–4. pmid:7877684
  2. 2. Desai D, Gu Y, Morgan DO. Activation of human cyclin-dependent kinases in vitro. Mol Biol Cell 1992;3:571–82. pmid:1535244
  3. 3. Connell-Crowley L, Solomon MJ, Wei N, Harper JW. Phosphorylation independent activation of human cyclin-dependent kinase 2 by cyclin A in vitro. Mol Biol Cell 1993;4:79–92. pmid:8443411
  4. 4. Bloom J, Cross FR. Multiple levels of cyclin specificity in cell-cycle control. Nature Reviews Molecular Cell Biology 2007;8:149–60. pmid:17245415
  5. 5. Gould KL, Moreno S, Owen DJ, Sazer S, Nurse P. Phosphorylation at Thr167 is required for Schizosaccharomyces pombe p34cdc2 function. EMBO J 1991;10:3297–309. pmid:1655416
  6. 6. Krek W, Nigg EA. Cell cycle regulation of vertebrate p34cdc2 activity: identification of Thr161 as an essential in vivo phosphorylation site. New Biol 1992;4:323–9. pmid:1622929
  7. 7. Gu Y, Rosenblatt J, Morgan DO. Cell cycle regulation of CDK2 activity by phosphorylation of Thr160 and Tyr15. EMBO J 1992;11:3995–4005. pmid:1396589
  8. 8. Solomon MJ. Activation of the various cyclin/cdc2 protein kinases. Current Opinion in Cell Biology 1993;5:180–6. pmid:8507489
  9. 9. Ducommun B, Brambilla P, Félix MA, Franza BR, Karsenti E, Draetta G. cdc2 phosphorylation is required for its interaction with cyclin. EMBO J 1991;10:3311–9. pmid:1833185
  10. 10. Stevenson LM, Deal MS, Hagopian JC, Lew J. Activation Mechanism of CDK2: Role of Cyclin Binding versus Phosphorylation. Biochemistry 2002;41:8528–34. pmid:12081504
  11. 11. Lapenna S, Giordano A. Cell cycle kinases as therapeutic targets for cancer. Nat Rev Drug Discov 2009;8:547–66. pmid:19568282
  12. 12. Peyressatre M, Prével C, Pellerano M, Morris MC. Targeting cyclin-dependent kinases in human cancers: from small molecules to Peptide inhibitors. Cancers (Basel) 2015;7:179–237. pmid:25625291
  13. 13. Asghar U, Witkiewicz AK, Turner NC, Knudsen ES. The history and future of targeting cyclin-dependent kinases in cancer therapy. Nat Rev Drug Discov 2015;14:130–46. pmid:25633797
  14. 14. Jeffrey PD, Russo AA, Polyak K, Gibbs E, Hurwitz J, Massagué J, et al. Mechanism of CDK activation revealed by the structure of a cyclinA-CDK2 complex. Nature 1995;376:313. pmid:7630397
  15. 15. Russo AA, Jeffrey PD, Pavletich NP. Structural basis of cyclin-dependent kinase activation by phosphorylation. Nature Structural Biology 1996;3:696. pmid:8756328
  16. 16. De Bondt HL, Rosenblatt J, Jancarik J, Jones HD, Morgan DO, Kim SH. Crystal structure of cyclin-dependent kinase 2. Nature 1993;363:595–602. pmid:8510751
  17. 17. Brown NR, Noble ME, Endicott JA, Johnson LN. The structural basis for specificity of substrate and recruitment peptides for cyclin-dependent kinases. Nat Cell Biol 1999;1:438–43. pmid:10559988
  18. 18. Brown NR, Noble ME, Lawrie AM, Morris MC, Tunnah P, Divita G, et al. Effects of phosphorylation of threonine 160 on cyclin-dependent kinase 2 structure and activity. J Biol Chem 1999;274:8746–56. pmid:10085115
  19. 19. Adams JA. Kinetic and Catalytic Mechanisms of Protein Kinases. Chem Rev 2001;101:2271–90. pmid:11749373
  20. 20. Lassila JK, Zalatan JG, Herschlag D. Biological phosphoryl-transfer reactions: understanding mechanism and catalysis. Annu Rev Biochem 2011;80:669–702. pmid:21513457
  21. 21. Bao ZQ, Jacobsen DM, Young MA. Briefly bound to activate: transient binding of a second catalytic magnesium activates the structure and dynamics of CDK2 kinase for catalysis. Structure 2011;19:675–90. pmid:21565702
  22. 22. Jacobsen DM, Bao Z-Q, O’Brien P, Brooks CL, Young MA. Price To Be Paid for Two-Metal Catalysis: Magnesium Ions That Accelerate Chemistry Unavoidably Limit Product Release from a Protein Kinase. J Am Chem Soc 2012;134:15357–70. pmid:22891849
  23. 23. Cheng Y, Zhang Y, McCammon JA. How Does the cAMP-Dependent Protein Kinase Catalyze the Phosphorylation Reaction: An ab Initio QM/MM Study. J Am Chem Soc 2005;127:1553–62. pmid:15686389
  24. 24. De Vivo M, Cavalli A, Carloni P, Recanatini M. Computational Study of the Phosphoryl Transfer Catalyzed by a Cyclin-Dependent Kinase. Chemistry—A European Journal 2007;13:8437–44. pmid:17636466
  25. 25. Smith GK, Ke Z, Guo H, Hengge AC. Insights into the Phosphoryl Transfer Mechanism of Cyclin-dependent Protein Kinases from Ab Initio QM/MM Free-energy Studies. J Phys Chem B 2011;115:13713–22. pmid:21999515
  26. 26. De Vivo M, Ensing B, Dal Peraro M, Gomez GA, Christianson DW, Klein ML. Proton shuttles and phosphatase activity in soluble epoxide hydrolase. J Am Chem Soc 2007;129:387–94. pmid:17212419
  27. 27. Cavalli A, Vivo MD, Recanatini M. Density functional study of the enzymatic reaction catalyzed by a cyclin-dependent kinase. Chem Commun 2003;0:1308–9.
  28. 28. Dubecký M, Walter NG, Šponer J, Otyepka M, Banáš P. Chemical feasibility of the general acid/base mechanism of glmS ribozyme self-cleavage. Biopolymers 2015;103:550–62. pmid:25858644
  29. 29. Maestro (2014) Version 9.7, Schrödinger, LLC, New York, NY.
    • 30. Sastry GM, Adzhigirey M, Day T, Annabhimoju R, Sherman W. Protein and ligand preparation: parameters, protocols, and influence on virtual screening enrichments. J Comput Aided Mol Des 2013;27:221–34. pmid:23579614
    • 31. Schrödinger Suite 2014–1 Protein Preparation Wizard.
      • 32. Olsson MHM, Søndergaard CR, Rostkowski M, Jensen JH. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions. J Chem Theory Comput 2011;7:525–37. pmid:26596171
      • 33. Banks JL, Beard HS, Cao Y, Cho AE, Damm W, Farid R, et al. Integrated Modeling Program, Applied Chemical Theory (IMPACT). J Comput Chem 2005;26:1752–80. pmid:16211539
      • 34. Berendsen HJC, Postma JPM, van Gunsteren WF, Hermans J. Interaction Models for Water in Relation to Protein Hydration. In: Pullman B, editor. Intermolecular Forces: Proceedings of the Fourteenth Jerusalem Symposium on Quantum Chemistry and Biochemistry Held in Jerusalem, Israel, April 13–16, 1981, Dordrecht: Springer Netherlands; 1981, p. 331–42.
        • 35. Desmond Molecular Dynamics System (2013) Version 3.4, D. E. Shaw Research, New York, NY.
          • 36. Essmann U, Perera L, Berkowitz ML, Darden T, Lee H, Pedersen LG. A smooth particle mesh Ewald method. J Chem Phys 1995;103:8577–93.
          • 37. Martyna GJ, Tobias DJ, Klein ML. Constant pressure molecular dynamics algorithms. J Chem Phys 1994;101:4177–89.
          • 38. Nosé S. A unified formulation of the constant temperature molecular dynamics methods. J Chem Phys 1984;81:511–9.
          • 39. Tuckerman M, Berne BJ, Martyna GJ. Reversible multiple time scale molecular dynamics. J Chem Phys 1992;97:1990–2001.
          • 40. Martí S, Moliner V, Tuñón I. Improving the QM/MM Description of Chemical Processes: A Dual Level Strategy To Explore the Potential Energy Surface in Very Large Systems. J Chem Theory Comput 2005;1:1008–16. pmid:26641916
          • 41. Becke AD. Density-functional thermochemistry. III. The role of exact exchange. J Chem Phys 1993;98:5648–52.
          • 42. Jorgensen WL, Tirado-Rives J. The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin. J Am Chem Soc 1988;110:1657–66. pmid:27557051
          • 43. Field MJ, Albe M, Bret C, Martin FP-D, Thomas A. The dynamo library for molecular simulations using hybrid quantum mechanical and molecular mechanical potentials. Journal of Computational Chemistry 2000;21:1088–100.
          • 44. Foster JP, Weinhold F. Natural hybrid orbitals. J Am Chem Soc 1980;102:7211–8.
          • 45. Reed AE, Curtiss LA, Weinhold F. Intermolecular interactions from a natural bond orbital, donor-acceptor viewpoint. Chem Rev 1988;88:899–926.
          • 46. Gaussian 09, Revision C.01, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, 25 S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, Gaussian, Inc., Wallingford CT, 2009.
            • 47. Pérez-Gallegos A, Garcia-Viloca M, González-Lafont À, Lluch JM. SP20 Phosphorylation Reaction Catalyzed by Protein Kinase A: QM/MM Calculations Based on Recently Determined Crystallographic Structures. ACS Catal 2015;5:4897–912.
            • 48. Montenegro M, Garcia-Viloca M, Lluch JM, González-Lafont À. A QM/MM study of the phosphoryl transfer to the Kemptide substrate catalyzed by protein kinase A. The effect of the phosphorylation state of the protein on the mechanism. Phys Chem Chem Phys 2010;13:530–9. pmid:21052604
            • 49. Pérez-Gallegos A, Garcia-Viloca M, González-Lafont À, Lluch JM. A QM/MM study of the associative mechanism for the phosphorylation reaction catalyzed by protein kinase A and its D166A mutant. J Comput Aided Mol Des 2014;28:1077–91. pmid:25129483
            • 50. Díaz N, Field MJ. Insights into the phosphoryl-transfer mechanism of cAMP-dependent protein kinase from quantum chemical calculations and molecular dynamics simulations. J Am Chem Soc 2004;126:529–42. pmid:14719950
            • 51. Zhao Y, Truhlar DG. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor Chem Account 2008;120:215–41.
            • 52. Mildvan AS. Mechanisms of signaling and related enzymes. Proteins 1997;29:401–16. pmid:9408938
            • 53. Zhang J, Li C, Shi T, Chen K, Shen X, Jiang H. Lys169 of Human Glucokinase Is a Determinant for Glucose Phosphorylation: Implication for the Atomic Mechanism of Glucokinase Catalysis. PLOS ONE 2009;4:e6304. pmid:19617908
            • 54. Bevilacqua PC, Brown TS, Nakano S, Yajima R. Catalytic roles for proton transfer and protonation in ribozymes. Biopolymers 2004;73:90–109. pmid:14691943
            • 55. Perreault DM, Anslyn EV. Unifying the Current Data on the Mechanism of Cleavage–Transesterification of RNA. Angewandte Chemie International Edition in English 1997;36:432–50.
            • 56. Lopez X, Schaefer M, Dejaegere A, Karplus M. Theoretical Evaluation of pKa in Phosphoranes: Implications for Phosphate Ester Hydrolysis. J Am Chem Soc 2002;124:5010–8. pmid:11982365
            • 57. NBO version 3.1, Glendening ED, Reed AE, Carpenter JE, Weinhold F (1990).
              • 58. Vogt-Geisse S, Toro-Labbé A. The mechanism of the interstellar isomerization reaction HOC+→HCO+ catalyzed by H2: New Insights from the reaction electronic flux. J Chem Phys 2009;130:244308. pmid:19566154
              • 59. Duarte F, Vöhringer-Martinez E, Toro-Labbé A. Insights on the mechanism of proton transfer reactions in amino acids. Phys Chem Chem Phys 2011;13:7773–82. pmid:21442087
              • 60. Reed AE, Weinstock RB, Weinhold F. Natural population analysis. J Chem Phys 1985;83:735–46.
              • 61. Roston D, Demapan D, Cui Q. Leaving Group Ability Observably Affects Transition State Structure in a Single Enzyme Active Site. J Am Chem Soc 2016;138:7386–94. pmid:27186960
              • 62. Baxter NJ, Bowler MW, Alizadeh T, Cliff MJ, Hounslow AM, Wu B, et al. Atomic details of near-transition state conformers for enzyme phosphoryl transfer revealed by rather than by phosphoranes. PNAS 2010;107:4555–60. pmid:20164409
              • 63. Dewick Paul M. Essentials of Organic Chemistry: For Students of Pharmacy, Medicinal Chemistry and Biological Chemistry. John Wiley & Sons; 2013.
                • 64. Xie Y, Jiang Y, Ben-Amotz D. Detection of amino acid and peptide phosphate protonation using Raman spectroscopy. Anal Biochem 2005;343:223–30. pmid:16018962