The Mechanism of Phosphoryl Transfer Reaction and the Role of Active Site Residues on the Basis of Ribokinase-like Kinases

The role of ribokinase-like carbohydrate kinases consists in ATP dependent phosphorylation of small molecules containing hydroxymethyl group. Although they differ substantially in structural terms and exhibit a broad substrate specificity, some family-wide conserved features can be distinguished suggesting the common mode of action. 4-methyl-5--hydroxyethylthiazole kinase (Thz kinase) was chosen as a representative model and the mechanism proposed in X-ray crystal structure paper provided the basis for calculations. In particular, the possible role of several active site residues (Arg121 and Cys198 among others) and of the two magnesium ions was examined. Static and dynamic catalytic fields for the reaction were generated revealing the most favourable environment for the preferential transition state stabilization. An attempt to model the phosphoryl transfer reaction as well as to investigate the influence of the cysteine residue on the reaction course at the semiempirical PM3 level of theory was undertaken.


Introduction
Structure of Bacillus subtilis YXKO -hypothetical protein with unknown function obtained recently in Structural Genomics project [1], revealed structural homology to ribokinase-like family of carbohydrate kinases (RK-kinases) and raised some questions about catalytic mechanism.
Kinases are a ubiquitous group of enzymes that facilitates phosphoryl transfer reaction from International Journal of

Molecular Sciences
ISSN 1422-0067 © 2004 by MDPI www.mdpi.org/ijms/ a phosphate donor (usually ATP) to an acceptor substrate.The largest and the most extensively studied is the family of protein kinases, but other groups of kinases, neither so well characterized nor investigated, also exhibit very interesting features.
Despite of an observed diversity, the overall geometry of substrate binding pocket is similar and hydroxymethyl group of phosphate acceptor molecule is positioned in nearly identical location.Sequence alignment performed on the basis of structure comparison also reveals some structurally significant motifs that are highly conserved [4].Thus, it is of interest not only to understand how phosphoryl transfer occurs in this particular group of enzymes, but also to understand all the differences and similarities between separate RK-like kinases.Taking the fact that conserved residues play a vital role in substrates binding and performing catalysis into consideration, a common catalytic mechanism was suggested for the ribokinase family, which implicates both modes of kinase catalysis, that is substrate activation and transition state stabilization [3][4][5]7].
Generally, phosphoryl transfer proceeds as an ordered process without involvement of a phosphoenzyme intermediate.The hydroxyl group of acceptor molecule is activated by a basic group (the suggested catalytic base utilized in most cases in deprotonation of hydroxyl is the carboxylate of Asp residue) and then its oxygen atom performs a direct nucleophilic attack to the  phosphate group of ATP in an in-line displacement mechanism.Some key interactions are involved in the correct positioning and stabilization of the negatively charged phosphoryl group and resulting pentacovalent transition state including presence of oxyanion hole, utilizing magnesium ion coordination, interaction of the guanidinum moiety of Arg with  phosphate and the positive helix dipole.
Although some facts can support an associative S N 2 type mechanism, no computational analysis at the quantum mechanical level have been made for this type of kinases.There are some more not fully clarified aspects of catalysis including the lack of suggested catalytic base in case of 4-methyl-5-hydroxyethylthiazole kinase (Thz kinase, ThiK) and 2-methyl-4-amino-5-hydroxymethylpirimidine kinase (HMPP kinase [8]).Questions addressed above, together with the availability of some ThiK mutants analysis, resulted in the choice of Thz kinase from Bacillus subtilis [7] as a model for investigation of some aspects of RK-like kinases catalytic mechanism.
ThiK is involved in thiamin biosynthesis and carbohydrate metabolism.The reaction facilitated by this particular kinase -phosphoryl transfer from ATP to 4-methyl-5--hydroxyethylthiazole (Thz), is illustrated in Figure 1.The crystal structures of the enzyme and its complexes: enzyme/Thz and enzyme/ATP/Thzphosphate seem to confirm the mechanism suggested for RK-like kinases: an in-line nucleophilic attack of oxygen atom from hydroxyl group of thiazole on  phosphorus of ATP.However, in case of ThiK there is no family-wide conserved Asp residue which is believed to serve as alcohol activating base.Instead of Asp, ThiK has a conserved Cys residue at the same position (Cys198) and it was experimentally shown that mutation of this cysteine to aspartic acid results in 9-fold higher enzymatic activity [7].However, the substitution of Cys by serine and alanine does not implicate the loss of enzymatic activity (the specific activity of corresponding mutants is 20 and 40% that of the native protein).
Although ThiK, like the rest of RK-like family members, utilizes a magnesium ion and an oxyanion hole to stabilize the negatively charged phosphoryl group, there is an additional electrostatic interaction between  phosphate of ATP and guanidinium moiety of residue Arg121.Arginine (or equivalent lysine) at this position is absolutely conserved within ThiK kinase family members [7] and has been recently found also in YXKO -putative RK-like kinase from B. subtilis [1].
The role of these particular residues, as well as some other aspects of phosphoryl transfer reaction catalyzed by ThiK as a representative of RK-like kinases are the subject of this study.In particular, some computer aided attempts to model the enzymatic reaction pathway and reveal the optimal catalytic environment for the process were carried out.

Preparation of the Starting Geometries
Crystal structures of enzyme bound to Thz and ATP/Thz-phosphate were obtained from Protein Data Bank (access codes 1C3Q and 1ESQ, respectively).The first one (the wild-type enzyme complexed with one of its substrates -Thz) is determined at the higher resolution (2.0Å comparing to 2.5Å) and contains no missing protein chain residues, however ATP molecule, as well as magnesium ions, are not included.The other crystal structure (the Cys198 to Ser mutant) incorporates substrate (ATP) and product (thiazole phosphate) molecules simultaneously, together with two magnesium ions.Thus, the two superimposed crystal structures were merged to reconstruct the entire active site.Two subunits with one substrate binding pocket located at their interface were chosen for simulations.Initial model was constructed by removing the phosphate from thiazole moiety and rotating the  phosphate of ATP into the active site.The one of magnesium ions that bridges  and  phosphates oxygen atoms is also coordinated through two water molecules to the aspartate and glutamate side chains (Asp 94 and Glu 126, respectively), whilst the second one is found to be complexed only with ATP (namely  and  phosphates) and does not interact with protein residues atoms.In case of cAMPdependent protein kinase (protein kinase A, PKA) corresponding magnesium ion is reported to inhibit the reaction but also to increase the thermal stability of its catalytic subunit when measured in the presence of ATP or ADP [9].According to recently reported tyrosine kinase (another protein kinase) mechanism [10], only one of the two magnesium ions present in the active site is involved in phosphoryl transfer reaction.Due to unclear significance of discussed magnesium ion, two models: with a single and with both magnesium ions were considered in further calculations.Crystallographic water molecules in the proximity of an active site were also included.

Molecular Mechanics Optimization
Optimization of protein complexed with ligands was performed using the molecular modeling package Insight II [11].Hydrogen atoms were built in Biopolymer module (assuming neutral pH).Additionally, the protonation state of ligands and active site residues was verified and corrected.Magnesium ions were included as divalent cations and ATP molecule charge was set to -4.In aqueous solution under neutral pH the triphosphate tail exists with the charge of -4, especially when ATP is complexed with magnesium ions [12].The difference in the conditions between protein cavity and water environment increases the number of possible protonation states, in most calculations, however, only fully deprotonated phosphoryl groups are considered [13,14].This investigation is also limited to such a case.
Since it was reported [15,16] that not only the solvent excluding optimization but also the choice of solvation model (namely fully macroscopic versus fully atomistic modeling approaches) affects crucial electrostatic interactions resulting in improper protein fold, the entire model was first solvated in 5 Å thick layer of water molecules and then submitted to optimization within Discover module.The Class II consistent quantum mechanics-based CFF97 forcefield was the source of parameters for potential functions.The energy minimization procedure was carried out in several stages to allow gradual relaxation of the system.At each stage, the steepest descent method was followed by conjugate gradient minimization algorithm.First, the positions of hydrogen atoms were optimized with fixed coordinates of remaining heavy atoms.Next, all water molecules were allowed to move and, finally, the protein and ligands heavy atoms' coordinates were refined by progressive decrease of tethering force.The final maximum energy derivative was 0.01 kcal/(molÅ).
The structures of substrates including one magnesium ion were extracted from optimized model and served as input geometry for further optimization.At this stage of investigation the model with two magnesium ions has not yet been included.

Quantum Mechanics Calculations
A model reaction system for the quantum mechanics calculations was then simplified: the adenine and the ribose moieties, as well as the aromatic ring of Thz alcohole, were not included.Instead of, the terminated form of ATP, methyl triphosphate and ethanol were used.Resulted model was optimized in the gas phase using Gaussian03 [17] at semiempirical (PM3 [18]) and ab initio (Hartree-Fock/STO-3G) levels of theory.Searches for transition state by means of Synchronous Transit-Guided Quasi-Newton (STQN) method [19] and optimization of approximate models to transition state instead of local minimum were unsuccessful: the normal mode of vibration in resulted structures associated with the negative force constant was not responsible for the transfer of the phosphoryl group.For further simulations transition state proposed for similar reaction catalyzed by protein kinase A [20] was used.

Static and Dynamic Catalytic Fields
Differential transition state stabilization (DTSS) approach [21] was applied to reveal the optimal configuration of catalytic environment.The catalytic activity can be considered in terms of activation barrier change , which is understood as the difference between activation energies for the reaction in the presence and in the absence of catalyst (in the gas phase) [21,22].From mathematical point of view,  is equivalent to the difference in energies of interaction between catalytic environment and transition state with respect to interaction of the same molecular environment with the substrates (or products).Electrostatic effects are one of the most significant contribution to interaction energies (especially in case of polar or charged systems) [22].Thus, assuming the dominant role of electrostatics in environmental effects and employing electrostatic approximation for the reactants only, it is possible to estimate the catalytic or inhibitory effect of the enzyme active site environment.The static [21] and dynamic [23] characteristics of the most favourable catalytic surroundings presented as a differential maps of electrostatic molecular potentials and vector field may be simply expressed as Δ s =V TS −V SC and −grad Δ d =  E TS −  E SC (where TS and SC denotes transition state and substrates complex, respectively;  E designates electric field vector; V and  E are calculated as an expectation value of a corresponding operator action on the wavefunction).The sign of electrostatic potential and the direction of electric field on the van der Waals contour is then reverted in order to reflect the complementary influence of environment.
Analysis of static and dynamic catalytic fields on the van der Waals surface grid surrounding superimposed substrates (excluding magnesium ion) and a putative transition state were performed.Electrostatic potential and electric field were calculated using Gaussian03 at HF/ STO-3G.

The Reaction Pathway Mapping
Calculations of a reaction pathway along an internal reaction coordinate defined as the distance between  phosphorus and oxygen atom of hydroxyl group were carried out using MOPAC7 [24] together with subroutine DRIVER [25].The distance was gradually decreased by 0.05Å with each step and the model was optimized with exception of the reaction coordinate.The position of methyl triphosphate carbon atom was frozen.The energy (heat of formation) of the minimized geometry in each step was calculated.As in previous case PM3 Hamiltonian was applied and model system was limited to that described for quantum mechanical calculations.Program Triton 3.0 [26,27] was employed to analyse the results.To reveal the influence of Cys residue on the reaction course, the simulation was repeated in the presence of cysteine (the coordinates of backbone atoms were fixed, while the side chain atoms were allowed to move).

Geometry Optimization
The two energy minimized models obtained from MM calculations were superimposed in order to reveal the influence of allowing for the second magnesium ion in the initial structure (Figure 2).Although the overall geometry of Thz molecule is nearly the same (with the RMSD value of the order of 0.1Å, when calculated for all atoms of superimposed thiazole structures only), there are some changes in ATP molecules conformations, but they are applied almost exclusively to triphosphate tail (corresponding RMSD value is 1.1 Å).The relative position of active site residues is also slightly different, especially when concerning the charged or polar residues that make the closest contact with substrates (Asp94, Arg121, Glu126, Cys198).More detailed analysis of the magnesium ions proximity is presented in Figure 3.The values of the distances between the Mg 2+ and oxygen atoms of water (in the range of 2.0-2.1 Å in all cases) molecules suggest a presence of interactions between these atoms.Except two water molecules coordinated to the second magnesium ion (the one bridging  and  phosphates) of the two Mg 2+ containing model, the origin of all concerned water molecules is the crystallographic structure.The remaining two water molecules were included as a solvent and then moved in the Mg 2+ proximity in the course of optimization.There are also close contacts between hydrogen atoms of water molecules coordinated to the Mg 2+ bridging  and  phosphates and negatively charged side chains of Asp94 and Glu126 in both models.The distances between reacting phosphorus atom of  phosphate and the thiazole oxygen atom are 3.6 and 3.9 Å, when comparing the two models.
Further gas phase optimization of simplified model containing methyl triphosphate, ethanol, single Mg 2+ bridging  and  phosphates and cysteine residue did not reveal any significant changes.As previously, magnesium ion is 1.8Å distant from coordinated  and  phosphate oxygen atoms.The distance between  phosphorus atom and hydroxyl oxygen atom slightly increased from 3.6 to 3.8Å Resulting structure (Figure 4) was utilized as a starting point for following calculations.

Static and dynamic catalytic fields
The comparison of two-dimensional differential electrostatic potential and electric field vectors plots around superimposed molecules with the corresponding active site geometries is presented in Figure 5.The magnitude of signs and arrows is proportional to the differential values and reflects the activation energy barrier lowering induced by the corresponding unit charge or the positive probe charge movement coupled with the reaction coordinate, respectively.The larger is the sign or arrow magnitude, the bigger is the coupled activation barrier lowering.
It is remarkable, that the positions of both magnesium ions as well as most of the charged amino acid residues in the active site (Asn25 and Arg121) coincide with the predicted topology of the static catalytic fields with exception of Asp94 and Glu126 screened by one of magnesium ions.

The Reaction Pathway Mapping -Preliminary Results
MOPAC/DRIVER [24,25] mapping of the reaction pathway results for two model systems (with and without cysteine residue) are compared in Figure 6.Despite of decreasing the distance corresponding to the reaction coordinate and -as a result -creating the bond between reacting atoms, the  phosphate remained bound to the rest of ATP fragment.However, the influence of the Cys198 presence is apparent and implicates the lowering of the reaction energy barrier.Also the energy level corresponding to the final structure with cysteine is located beneath the one that corresponds to the model where cysteine is not included.Besides, the maximum and minimum in energy in case of cysteine containing model are reached earlier and the value of the reaction coordinate for the minimum energy structure (~2.1Å) is too large to indicate the final product forming.Moreover, it is close to the value of the reaction coordinate distance in maximum energy structure of a model without cysteine (namely ~2.2Å).Even though the minimum energy structure of this model was obtained for ~1.8Å, it is still more than expected value.In both cases proton transfer from hydroxyl group to the magnesium coordinated oxygen atom occurs.It should be noted, that the geometries of models corresponding to the energy maxima were not verified in terms of vibrational frequencies, thus they are considered as the putative transition state structures only.Figure 6.The energy (heat of formation) changes along the reaction coordinate (distance of  phosphorus atom from the nucleophilic oxygen atom of hydroxyl group); broken and full horizontal lines correspond to the substrates only and substrates together with cysteine containing models, respectively.Geometries of investigated models for the marked levels of energy are also presented.In parentheses are given the related reaction coordinate values.For both models the initial reaction coordinate value is 4Å.

Discussion and Concluding Remarks
Molecular mechanics optimization of starting geometries resulted in reasonable structures with a good correspondence to an active site model proposed in [7].Especially, the observed water molecules network around magnesium ions that interacts also with two negatively charged residues (in case of Mg 2+ included in both models) is worth noticing.However, since magnesium ion is generally hexa coordinated, the solvation shell of the magnesium ion bridging  and  phosphates may still be uncomplete (only five ligands are observed: two oxygen atoms of phosphates and three water molecules).An existence of a set of conserved water molecules was reported, for example, in case of protein kinase A [28], where comparison of known crystal structures revealed the presence of six structured water molecules interacting with ATP, two magnesium ions (each of them coordinated with six ligands) and such an active site residues as aspartate, glutamate or lysine.Although the optimized distances between water molecules oxygen atoms and magnesium ions are similar to that proposed for protein kinase A (of the order of 2.1Å), the contacts made by Mg 2+ and coordinated phosphate oxygen atoms (1.8Å) are too close (expected value is in the range of 2.0 -2.2Å [12]).
Generation of static and dynamic catalytic fields revealed the general characteristics of the electrostatic environment, which is favourable for catalysis by stabilization of the transition state rather than the substrates.Comparison of the optimal catalytic environment obtained only from the reactants structures with the actual geometry of ThiK active site indicates that catalytic fields are generally in a good agreement with the positions of the most critical moieties in the substrate binding cavity.The location of positive charges is most clearly seen.According to this, both magnesium ions seem to be important for preferential transition state stabilization and activation barrier lowering.This implicates a necessity of taking the second Mg 2+ into account in further calculations.The presence of Arg121 is also confirmed as being important not only to compensate the negatively charged  phosphate of substrate ATP but also as a possible charge stabilizer in case of transition state.
There is no negatively charged residue in the proximity of attacking hydroxyl group in wild-type enzyme, however, proposed catalytic environment suggesting the present of a negative charge in this location could provide an explanation for the higher specific activity of C198D mutant.One way to verify such a hypothesis could be performing in silico mutations and comparing the results of reaction pathway mapping for a series of mutants.
As a conclusion, monitoring of the electrostatic properties changes during the reaction provides a simple but effective method to elucidate the essential electrostatic interactions between an enzyme and a substrate and can assist in the identification of the catalytically important active site residues, especially when considering charged or polar moieties.Moreover, such an approach is not limited to only those cases, when an accurate transition state for a given reaction is available -utilization of an approximate structure can also result in reasonable estimation of catalytic fields.
Preliminary investigation of reaction pathway seems to confirm the influence of Cys198 residue on the phosphoryl transfer process: the energy barrier between putative transition states (assigned as structures corresponding to maximum of energy changes curve) and starting geometries is lowered at approximately 9 kcal/mol in the presence of cysteine.From the other side, since there is neither bond creating nor proton/charge transfer between Thz and Cys198, its role may consist largely in proper alignment of attacking nucleophilic group (as suggested in [7]).Performed simulation resembles suggested associative type of phosphoryl transfer mechanism including pentacovalent geometry of  phosphorus atom.However, obtained results should be considered as a starting point for further study and a guideline for investigation to be carried out with larger models and more sophisticated methods, for example ab initio QM/MM.

Figure 2 .
Figure 2. Superimposition of MM optimized ThiK active site fragments (substrates and the side chains of the most important amino acid residues).The model containing one magnesium ion is presented in stick representation, white colour.Only heavy atoms of all molecules are shown.

Figure 3 .
Figure 3.The closest contacts made by magnesium ions (presented as green balls).For distinction, the two non-crystalline water molecule are shown in stick representation.The missing values of marked distances are described in the text.

Figure 4 .
Figure 4. Gas phase optimized geometry of investigated models (with cysteine residue shown in stick representation).

Figure 5 .
Figure 5. Symbolic representations of the differential electrostatic properties for two perpendicular planes (the top of the figure).Superimposed structures of both substrates (dotted lines, position of bonds only) and transition state (full lines) are presented.The bottom part of the figure: substrate binding pockets presented from views corresponding to plots (amino acid residues presented in stick representation).