Exploring the Mechanism of Shikimate Kinase through Quantum Mechanical and Molecular Mechanical (QM/MM) Methods

Abstract

The chemical step of Shikimate Kinase Helicobacter pylori, involving the transfer of a phosphoryl group, has been studied by using quantum mechanical and molecular mechanical (QM/MM) methods. Understanding the mechanism of this chemical step, present in bacteria and other microorganisms but absent in humans, can lead to the development of novel drugs for the treatment of common diseases caused by those pathogenic organisms. Different mechanisms including associative, dissociative, and concerted have been proposed up to now but there is not a consensus on the type of pathway that the reaction follows. Herein, we found that the mechanism has features from the associative and concerted types. An analysis of the free energy landscape of the chemical step reveals that the reaction is a two-step process without a well-defined intermediate state.

1. Introduction

Shikimate kinase (SK) catalyses the conversion of ATP and shikimic acid (SKM) to ADP and shikimate-3-phosphate (S3P) by transferring the terminal phosphoryl group of ATP to SKM. This phosphorylation reaction is an essential step in the overall shikimate pathway whose purpose is to produce the chorismate intermediate, a compound that is further processed in a chain of enzymatic reactions, that results in the synthesis of aromatic residues. Because this catalytic reaction is only found in microorganisms such as bacteria [1], for instance, Helicobacter pylori, which is responsible for gastrointestinal ulcers [2], it makes this reaction pathway a good target for novel drugs to inactivate such organisms. Several crystal structures of SK have been published up to now in the reactant state (RS) or product state (PS) [3,4] but a structure containing a transition state (TS) or TS mimic, which is crucial for drug development, has not been reported yet.

Depending on the surroundings, the phosphorylation reaction can follow any continuous pathway which is bounded by two extreme cases [5,6], also known as associative and dissociative mechanisms. In the associative mechanism, the bond between the leaving oxygen atom and the leaving phosphorus atom (P-O_lg) cleaves after the nucleophilic attack has started and the new bond between P and the nucleophilic group is formed (P-O_nuc). In the dissociative mechanism, the P-O_lg bond of the leaving group cleaves before the new bond P-O_nuc has been formed, see Figure 1 for an explanation of these distances. The associative mechanism is then characterised by a pentavalent phosphorene intermediate while the dissociative mechanism is by a planar metaphosphate intermediate. In between these two mechanisms, one can find a special case called the concerted mechanism where bond cleaving and forming occur simultaneously and the associative/dissociative character depends on the bond orders in the transition state. The reaction mechanism of SK has been studied previously both experimentally [2] and theoretically [7], however, there is not a consensus about the actual pathway for this reaction. For instance, based on the long distances between reactants observed in molecular dynamics simulations, Prado et al. [2] suggested that the reaction involves a dissociative mechanism. This finding was also supported by the fact that the catalytic turnover was achieved in experiments with ADP/metaphosphoric acid (an ATP-mimetic) in the absence of the native reactant ATP. This showed that a tight bond between the leaving group and ADP was not necessary. Contrary to the previous result, Hartmann et al. [3] proposed previously an associative reaction mechanism based on a geometric distance analysis from a model structure obtained from a series of crystal structures along the catalytic step. Recently, Yao et al. [7] proposed a concerted one-step reaction with a loose TS, by using QM/MM simulations. This was the first time that a model for the TS was proposed.

The difference in the previous proposed results regarding the mechanism of the SK’s chemical step motivated us to further investigate it by using different levels of theory within the QM/MM [8,9,10] framework. They allowed us to compute free energies along a predefined reaction coordinate which in turn helped us to assess the type of mechanism involved with particular emphasis on the TS.

2. Materials and Methods

Helicobacter pylori SK in complex with shikimate-3-phosphate (S3P) and adenosine diphosphate (ADP), PDB ID 3MUF [4], was used as the initial structure of SK, S3P and ADP. One Mg²⁺ atom, present in other protein kinases [11], was manually added between ADP and S3P. In addition to this, four coordinating water molecules were added around the Mg²⁺ atom to achieve its six-fold coordination. CHARMM package (versions 40a2 and 43b1) [12] was used to build the initial system which contained the protein, ADP, S3P, Mg²⁺ ion, 59 crystal waters, 16 Na⁺ ions, 13 Cl⁻ ions, and 5184 bulk water molecules resulting in a system of 18,447 total number of atoms. The system was simulated as a rhombic dodecahedron (RHDO) box with a length of 63.6 Å. The CHARMM27 [13,14,15] force field was employed to simulate the protein chain, waters and ions. Water molecules were simulated with the TIP3P [16] model. SHAKE [17] constraints were applied to all bonds that included hydrogens. Electrostatic interactions were solved with the 6th-order particle mesh Ewald (PME) summation technique [18] with a 10 Å real space cut-off. The van der Waals interactions were truncated with a switching function between 9 Å and 10 Å.

The MM part in our QM/MM simulations was performed with the CHARMM engine while the QM part was solved with either the MNDO97 [19,20] CHARMM module or the Gaussian 09 [21] package through the CHARMM interface [10]. The QM region included ADP, S3P, four coordinating waters, Mg²⁺, and the terminal atoms of the following residues: K14 (cut at C^δ atom), protonated D33 (cut at C^α), R57 (cut at C^γ atom), R107 (cut at C^γ atom), and R116 (cut at C^γ atom). The cutoff scheme for treating QM boundary atoms was the generalized hybrid orbital (GHO) method [22]. Although residues R57 and R107 are not close to the chemically important region, we included them in the QM region to zero out the total charge of it. The total number of QM atoms was 131. In the MNDO simulations, the Mg²⁺, S3P and ADP molecules were modelled with the AM1/d-PhoT [23] set of parameters. As for the Gaussian 09 simulations, the level of theory was B3LYP [24,25] with the 6-31+g(d) basis set. The extended Lagrangian formulation with dissipation in the Born-Oppenheimer approximation [26,27] was used in the MNDO simulations to speed up the convergence of the self-consistent procedure.

The initial structure (in the PS) was first minimised using the ABNR method [28] during 1000 steps while keeping the protein chain, ADP, S3P, and Mg²⁺ ion fixed (excluding all hydrogens). A subsequent minimisation of 500 steps was done where the constraints on the protein chain were relaxed. The generation of the initial pathway started with the latter structure and subsequent structures were obtained from the previous one by gradually reducing the value of the reaction coordinate (RC) which was computed as:

$$\mu =r_{P-O_{lg}}-r_{P-O_{nuc}} (\AA )$$

(1)

And whose initial value was µ = 2.5 Å. Each reduction of the RC was followed by a minimisation procedure of 5 and 500 steps for B3LYP and MNDO methods, respectively. Two pathways were generated by using B3LYP (500 windows) and the computationally cheaper semiempirical MNDO methods (50 windows). Free energy calculations were performed on the pathway obtained with the MNDO method by using the umbrella sampling [29] (US) technique. These finite temperature simulations were performed at 298.15 K in the constant pressure (CPT) ensemble with the Nosé–Hoover Langevin piston method [30,31].

The dissociative character of the TS, calculated as DC = (1 − n) × 100%, was derived from Pauling’s Equation [32,33,34]:

$$D\left(n\right)=D\left(1\right)-0.6\log \left(n\right)$$

(2)

where D(1) is the single bond distance for a P-O bond (1.73 Å) and D(n) is the average of the distances at the TS.

3. Results

3.1. Details of the Initial Structure in the PS

The binding site of the initially minimised structure with the MNDO method, see Figure 1, shows that the products ADP and S3P have a network of hydrogen bonds with the nearby residues K14 and R116 and therefore they function as a bridge between the products. This network seems to be essential for the catalytic efficiency of the enzyme as demonstrated by the dramatic reduction of activity upon mutations of these residues [4,35]. In this structure, D33 is protonated at the PS and forms a hydrogen bond with the nucleophilic oxygen (O_nuc). The Mg²⁺ ion helps to stabilise the negatively charged leaving group and to keep the leaving group in place.

3.2. Reaction Mechanism of the Chemical Step

The B3LYP-QM/MM minimised pathway shows that the distances between P-O_lg and P-O_nuc gradually change as was found by Yao et al. [7] with the crossing point of the curves occurring at µ = 0 Å. The distances for P-O_lg and P-O_nuc at this point were both equal to 2.64 Å, while the ones reported in this reference had values of 2.37 Å and 2.57 Å, respectively, see Figure 2a. The application of Pauling’s formula to these sets of distances results in a dissociative character DC obtained from Equation (2) of 96.9%, in our case and a close result of 94.1% in Yao’s DFT calculations. We noticed that the dissociative character is in fact closer to our value in the case of Yao’s MP2 calculation because the bonds at the TS are slightly looser than in the DFT results cited above. We also noticed that although P-O_lg and P-O_nuc gradually change, the slope of the curves changes producing a bending feature near the crossing point.

This bending feature is better captured in a More O’Ferrall Jencks plot [36,37] as it can be observed in Figure 2b. A plot for the B3LYP (red line) result shows that starting from the RS the reaction follows an associative mechanism, but close to the crossing point µ = 0 Å the curvature of the P-O_lg vs. P-O_nuc curves changes. A similar pattern was observed with MNDO (black line) level of theory. This bending suggests that the mechanism is shifted towards a more concerted one depicted by the green dashed line in the same figure. It is interesting to notice that in the case of the MNDO result, the reaction pathway is even closer to the concerted mechanism than in the B3LYP case which is a consequence of the looser bonds obtained with the former method, for instance at the crossing point the P-O_lg and P-O_nuc distances were both equal to 2.64 Å and 2.9 Å for B3LYP and MNDO methods, respectively.

Regarding the crossing point, we performed a frontier molecular orbital analysis on the structure obtained with B3LYP method at µ = 0 Å, this can be seen in Figure 3. We found that the highest occupied molecular orbital (HOMO), Figure 3a, is localised in the β phosphate group of ADP but the expected location was at the level of O_nuc as this region is supposed to donate the highest-energy electrons for attacking the leaving group. This is not the case in the present reaction mechanism because the proton of SKM is transferred to D33 beyond the crossing point (µ > 0.42 Å) and the attacking oxygen O_nuc is not yet enabled, this is shown in Figure 2a (green line). The lowest unoccupied molecular orbital (LUMO), Figure 3b, is localised around the Mg²⁺ ion showing that this ion assists in the stabilisation of the negative charges (red blob) formed when reactants approach each other. All these findings suggest that the crossing point is not a real transition state in our case but only one intermediate in an ensemble of states which can be accessed due to the loose nature of the bonds. The latter can be assessed by turning on a thermostat (finite temperature) which will make the relevant bonds fluctuate and therefore prove the stability of the state at the crossing point.

3.3. Free Energy Landscape of the Chemical Step

In order to explore the behaviour of the proposed reaction mechanism, with particular emphasis on the crossing point, we performed finite temperature simulations and computed the potential of mean force (PMF) along the reaction coordinate defined in Equation (1), the results are displayed in Figure 4. The PMF shows that the reaction proceeds through two transition states TS1 (µ = −0.75 Å) and TS2 (µ = 0.44 Å) without a clear stable intermediate state between them. This finding discards the fully associative/dissociative mechanisms where a stable phosphorene/methaphosphate intermediate is required. The role of TS1 is to decrease the overall barrier height by functioning as an intermediate state before the reactants climb to TS2, which is the rate-limiting step with a barrier height of 10.3 kcal/mol. This overall barrier height is lower than the experimental value of 15.3 kcal/mol derived from the k_cat value reported previously [35] by using Eyring’s formula of transition state theory [38]. It is known that semiempirical methods tend to underestimate the energy barriers [7] which are usually corrected by applying a more accurate but also more computationally expensive level of theory. However, we have not attempted to correct energies in this work because our main interest is on the pathway rather than on the exact energies. The PMF together with the More O’Ferrall plot confirms the idea that the reaction mechanism shares some features of the associative and concerted mechanisms.

Representative structures of the transition state ensemble TS1 and TS2 are displayed in Figure 5a,b, respectively. This figure shows that in the TS1, the leaving group is closer to ADP than to SKM which is the opposite in TS2. The Mg²⁺ ion helps to keep in place the leaving group. In both TS1 and TS2, the proton at the nucleophilic oxygen has not been transferred yet which contributes to the bond elongations.

4. Discussion

The analysis of the P-O_lg and P-O_nuc distances along the reaction coordinates reveals that the reaction starts by following an associative mechanism. However, as it approaches the crossing point µ = 0 Å, the bond lengths become looser and the reaction tends to be shifted towards a more concerted mechanism, see Figure 2b. The crossing point showed a 96.9% dissociative character suggesting that the bonds at this point are loose. Thus, we claim that the reaction is neither fully associative nor concerted but that it has features from both types of mechanisms. These results can explain the findings from Prado et al. [2], that the catalytic turnover can proceed with the mimetic metaphosphoric acid, where a tight bond between the leaving phosphoryl group and ADP is not necessary. Our model can explain this finding because near the crossing point the bonds are loose and dissociative in nature. The previous finding of Hartmann [3], regarding the associative character of the reaction can be explained by our model because right at the beginning of the chemical step, the reaction follows an associative mechanism. Finally, the results of Yao [7], are also explained by our results as the reaction becomes concerted near the crossing point. Finite temperature simulations showed that the crossing point is not the TS of the reaction, but that it is separated into two TSs instead of the well-defined TS found in the previous reference.

The proposed two-step mechanism for the SK-catalysed phosphoryl transfer reaction is presented in Figure 6a–d. Reactants are kept in place by the two long arm residues K14 and R116 and the Mg²⁺ ion. Upon reaching the two TSs, K14 and R116 show high flexibility which is reflected in the displacement of their terminal groups. The flexibility allows these residues to form and break hydrogen bonds with the reactants with a low energy cost. The Mg²⁺ ion also assists in the rearrangement of the hydrogen bond pattern by changing its coordination, with the leaving and nucleophilic groups’ oxygen atoms, depending upon the reaction state. In Figure 6d we also see that D33 is the acceptor of the proton released by SKM.

5. Conclusions

In the present work, we analysed the SK-catalysed phosphoryl transfer reaction through QM/MM simulations at different levels of QM theory. We found that the reaction mechanism lies between the fully associative and dissociative mechanisms. Analysis of distances through the More O’Ferrall Jencks plot and the dissociative character showed that the reaction turns into a concerted mechanism near the crossing point of the P-O_lg and P-O_nuc distance curves. We also found that the reaction proceeds through two TSs. Our work reconciles in a single model previous results that suggested three different possible mechanisms.