Introduction

Understanding the underlying physical mechanism of the impressive rate accelerations achieved by enzymes has been a central question in biology. Transition-state theory has provided the fundamental framework since the rate of a reaction is dictated by the free-energy difference between the ground and transition states (Truhlar, 2015). A wealth of research has led to the general notion that enzymatic rate acceleration is due to the enzymes’ much higher affinity to the transition state (TS) relative to its substrates. This idea has been supported experimentally by the high affinities measured for transition-state analogues (TSA) relative to substrates, leading to the design of TSA as high-affine enzyme inhibitors (Lienhard, 1973; Schramm, 2007). Since the TS represents a maximum in the free-energy landscape, its structure cannot be determined directly due to its low probability and ephemeral nature, TSA-bound species have been extensively used as structural proxies.

The classic transition-state theory and structural visualization of enzyme/TSA complexes lead to a current view of quite unique structures on the dividing surface at the maximum in the free-energy landscape. In contrast, it is now generally accepted that proteins are large ensembles of conformations, a concept rooted in pioneering work by Frauenfelder and coworkers for protein function, in analogy to the ensemble/funnel concept of protein folding (Frauenfelder & Wolynes, 1985). Consequently, enzyme-substrate (ES) complexes - experimentally accessible as minima in the free-energy landscape - are composed of a vast ensemble of molecular configurations due to their macromolecular nature. Hence, catalysis by pathing through exclusive TS’s would pose a huge entropic bottleneck/energy barrier for enzyme-catalyzed reactions.

Motivated by this consideration, together with the marginal outcomes in current enzyme design relative to the catalytic power of naturally evolved enzymes, we set out to investigate the atomistic nature of transition states of the chemical step in an enzymatic cycle and its implication for rate acceleration by combining quantum-mechanics/molecular-mechanics (QM/MM) calculations and experiments. We take advantage of the groundbreaking work on TS theory by QM simulation developed in the field (Kamerlin et al., 2013; Masgrau & Truhlar, 2015; Warshel & Bora, 2016; Zinovjev & Tu ñón, 2017), well summarized in several reviews (Cui, 2016; Senn & Thiel, 2009; van der Kamp & Mulholland, 2013). We chose a phosphoryl-transfer (P-transfer) reaction due to its impressive enzymatic rate acceleration, as the uncatalyzed reactions have extremely high energy barriers (Kerns et al., 2015; Lassila et al., 2011). The chemical rationale for these high barriers and reasons as to why P-transfer reactions are ubiquitous in living organisms and are involved in almost all essential biological processes is elegantly discussed in a fundamental paper by Westheimer in 1987: “Why nature chose phosphates” (Westheimer, 1987). The central role of enzyme-catalyzed P-transfer reactions in biology (i.e., genetic code, energy storage, signaling, and compartmentalization) has led to extensive discussions and controversies of how enzymes catalyze these vital reactions so efficiently (Allen & Dunaway-Mariano, 2016; Kamerlin et al., 2013; Kerns et al., 2015; Lassila et al., 2011; Pabis et al., 2016; Westheimer, 1987). Here, we combine QM/MM calculations with experimental temperature- and pH-dependent kinetic studies and X-ray crystallography of adenylate kinase (Adk) to shed light on the question of transition-state structures for the chemical step in its catalytic cycle. While QM/MM methods have been extensively applied to study enzymatic reactions, yielding accurate and detailed descriptions of the molecular events during catalysis (Acevedo & Jorgensen, 2010; Cheng et al., 2005; Dinner et al., 2001; Hahn et al., 2015; Hirvonen et al., 2020; Karplus & Kuriyan, 2005; Lai & Cui, 2020a; López-Canut et al., 2011; Mones et al., 2013; Palermo et al., 2015; Rosta et al., 2014; Roston et al., 2018; Roston & Cui, 2016; Schwartz & Schramm, 2009; Senn & Thiel, 2009; Turjanski et al., 2009), we test here the reaction free-energy profiles from simulations subsequently by designed experiments. We deliver a revised notion of transition-state stabilization by the enzyme through our discovery of a structurally broad transition-state ensemble (TSE).

Results

Adk is an extensively studied phosphotransferase that catalyzes the reversible, roughly isoenergetic conversion of two ADP molecules into ATP and AMP (Figure 1a,b), thereby maintaining the cellular concentrations of these nucleotides (Dzeja & Terzic, 2009). It is an essential enzyme found in every cell and organism. Extensive structural, dynamic and kinetic studies led to a comprehensive picture of the overall Adk reaction mechanism and its underlying energy landscape (Figure 1a) (Beckstein et al., 2009; Berry et al., 1994, 2006; Henzler-Wildman et al., 2007; Kerns et al., 2015; Müller et al., 1996; Müller & Schulz, 1992; Wolf-Watz et al., 2004). The conformational change of the AMP- and ATP-lids closing and opening is crucial in the enzymatic cycle: Lid-closing positions the two substrates and the active-site residues for efficient chemistry (i.e., P-transfer) and prohibits the alternative, energetically favorable reaction of phosphoryl hydrolysis (Kerns et al., 2015). Lid-opening, essential for product release, and not P-transfer is the rate-limiting step in the catalytic cycle and Mg²⁺ accelerates both steps (Kerns et al., 2015). Here we investigate the actual P-transfer step computationally, followed by experimental testing of our computational results. This chemical step would take about 7000 years without the enzyme (Stockbridge & Wolfenden, 2009) compared to > than 5,000 s⁻¹ with the enzyme (Kerns et al., 2015).

Figure 1:

QM/MM calculations of the P-transfer step in Adk

In order to shed light on the mechanism by which Adk from Aquifex aeolicus catalyzes the phosphoryl-transfer reaction by more than 12 orders of magnitude (Kerns et al., 2015), we performed QM/MM simulations starting with two ADP molecules and Mg²⁺ bound in the active site. The starting structures for the simulations were prepared using the X-ray structure of Adk in complex with Ap5A (P1,P5-Di(adenosine-5′) pentaphosphate) and coordinated to Zn²⁺ (2RGX (Henzler-Wildman et al., 2007)). ADP-ADP coordinates were built using Ap5A as a template and Zn²⁺ was replaced by Mg²⁺. For the QM/MM simulations, the QM region was defined as the diphosphate moiety of both ADP molecules, the Mg²⁺ ion, plus the four coordinating water molecules. The rest of the system was described at molecular mechanics level using the AMBER ff99sb force field (Hornak et al., 2006) and was solvated with TIP3P water molecules (Figure 1 – figure supplement 1).

The equilibrated starting structures agree well with the X-ray structures of the enzyme bound to Mg²⁺/ADP (4CF7 (Kerns et al., 2015) Figure 1 – supplement 2). Steered molecular dynamics simulations were run in both the forward (ADP/ADP to ATP/AMP) and reverse direction (ATP/AMP to ADP/ADP) with Mg²⁺ present in the active site (Figure 1 – supplement 3). Since it is unknown whether the fully charged or monoprotonated nucleotide (on one β−ADP oxygen) state is the more reactive configuration, we performed QM/MM simulations for both cases. Free-energy profiles of the P-transfer were determined using Multiple Steered Molecular Dynamics and Jarzynski’s Relationship (Crespo et al., 2005; Jarzynski, 1997; Ramírez et al., 2014) (Figure 1c, and Figures 1 – supplements 3-5). The results reveal a much smaller free energy of activation for the fully charged nucleotide state (Δ_fG^‡ of 13 ± 0.9 kcal/mol) relative to monoprotonated state (Δ_fG^‡ of 23 ± 0.9 kcal/mol) (Table 1); hence, we conclude that this is the most reactive enzyme configuration.

Table 1.

When the calculations were repeated in the absence of Mg²⁺ (with the nucleotide monoprotonated, since the fully-charged nucleotide state prohibited the reaction), a large increase in the free energy activation barrier was observed relative to the Mg²⁺ bound system (Figure 1c and Figure 1 – supplement 5, Δ_fG^‡ of 34 ± 1.6 kcal/mol), in agreement with expectations from experiments (Kerns et al., 2015). We note that only a lower limit for the overall acceleration by Mg²⁺ (>10⁵-fold) could be estimated from published results (Kerns et al., 2015), since the P-transfer was too fast to be measured experimentally in the presence of Mg²⁺.

During revision, we tested the accuracy and robustness of our results from our original QM/MM level of theory by employing higher level pure DFT (PBE) free energy calculations similar to the approaches used by Ganguly et al. (Ganguly et al., 2020) (Figure 1 – supplement 6). The results are in excellent agreement with the free energy profiles obtained before (Figures 1c,e and Figures 2a,c) thus further supporting the computational results.

Figure 2:

The Transition-State Ensemble (TSE) – transferring phosphoryl group delocalized

The nature of the transition state of enzyme-catalyzed P-transfer reaction with respect to its associative and dissociative character has been of significant interest and heated debate (Kamerlin et al., 2013; Kamerlin & Wilkie, 2007; Lassila et al., 2011; Roston & Cui, 2016). The definition as well as theoretical and experimental approaches to distinguish between them are rooted in elegant and fundamental work on nonenzymatic P-transfer reactions (Duarte et al., 2015; Hengge, 2002; Hou et al., 2012; Hou & Cui, 2012; Kamerlin et al., 2013; Kerns et al., 2015; Kirby & Nome, 2015; Stockbridge & Wolfenden, 2009). For Adk, a concerted mechanism is observed with no intermediates (Figures 1c and 2). Instead of using the classification of associative versus dissociative character (Lassila et al., 2011), we rather apply the clear definition by Cui and coworkers of a tight versus loose TS for concerted P-transfer reactions as they are unambiguously related to the bond-order in the transition state (Lai & Cui, 2020a, 2020b; Roston et al., 2018; Roston & Cui, 2016). For Adk, a contraction of the active site is observed while passing through the TSE, as captured by the decrease in the distance between acceptor and leaving group during the transition (Figures 1d,e). In the absence of Mg²⁺, the acceptor and donor need to be within 4 Å for the reaction to occur (Figure 1d) whereas in the presence of Mg²⁺ an acceptor-donor distance of 4.5 Å is sufficient (Figure 1e). The character of the TSE can best be seen from the widely used Moore-O’Ferrall-Jencks diagram (Jencks, 2002; O’Ferrall, 1970) that plots the distance from the transferring phosphate to the oxygen of the donor and acceptor, respectively (Figure 2). The TSE changes from tight/synchronous without Mg²⁺ (meaning high bond-order for both bonds between transferring phosphate to the leaving group and the attacking nucleophile) into loose with Mg²⁺ (lower bond-order for these two bonds; Figures 2 and 3) (Lai & Cui, 2020a, 2020b; Roston et al., 2018; Roston & Cui, 2016).

Figure 3:

Importantly, from these two-dimensional (Figure 2) and the one-dimensional (Figure 1c) free energy plots, we noticed a novel striking feature: a large ensemble of conformations in the TS region with vast differences in the position of the transferring phosphate but equal values in free energy for the enzyme with Mg²⁺ (light blue dots in Figure 2a, yellow area in Figure 1c). In other words, the enzyme seems to operate with a wide TSE. In contrast, without Mg²⁺ fewer conformations seem to comprise the TSE. To dive more into this unexpected computational result and better visualize the TSE, we show zoom-ins of representative TSE snapshots and their superposition (Figure 3 and Figure 3 – supplement 1). Notably, the TSE conformations are distributed along the reaction coordinate. The difference in position of the transferring phosphate along the reaction coordinate within the TSE (about 1 Å with Mg²⁺; Figures 1c, 2a, and 3b) implies that the TSE contains many highly asymmetric conformations, meaning that the transferring phosphate can be much closer to the leaving oxygen than the attacking oxygen and vice versa (Figure 3d). This asymmetry is logically tied with nonplanar configurations of the transferring phosphate (Figure 3d and Table 1).

Comparison of the TSE snapshots for the fully active, Mg²⁺-bound enzyme obtained from QM/MM with an X-ray structure of the same enzyme in complex with the transition-state analogue (ADP-Mg²⁺-AlF_4—AMP) (Kerns et al., 2015) serves as an initial experimental validation of our simulations (Figures 3e,f and Figure 3 – supplement 2). At the same time, this comparison highlights the power of the QM/MM simulations to investigate the catalytic mechanism, as the TSA-bound X-ray structure seems to imply quite a unique TS conformation, in sharp contrast to the broad TSE discovered in our simulations.

A detailed analysis of the pathways with and without Mg²⁺ reveals well-known features from other P-transfer enzymes, such as coordinated changes of the Mulliken charge populations coupled with proton transfers (Figure 3 – supplements 3-5). Mg²⁺ achieves its catalytic effect on the chemical step playing several roles. It enables the donor and acceptor phosphates to come close enough to react, thus achieving a reactive conformation. In the absence of Mg²⁺, a proton must bridge the terminal phosphates to prevent the strong electrostatic repulsion between the two nucleotides. In addition, Mg²⁺ stabilizes the charge of the transferring phosphate and thereby lowering the enthalpic barrier. Most interestingly, the presence of a cation seems to result in a wide TSE. Structures of the transition-state ensemble show that Mg²⁺ is optimally coordinated to the transferring phosphates and water molecules (Figure 3g).

Besides the multifaceted role of Mg²⁺, our simulations provide insights into the function of the fully-conserved arginine residues in the active site for lowering the activation barrier for this chemical reaction (Figure 3 – supplements 6-8). R85 appears to arrange the β-phosphates of the ADPs in the proper position for the reaction to be started. R150 and R161 are involved in anchoring ATP and AMP residues in the backward reaction. These last two residues and R85 interact with the Pα of AMP, stabilizing its negative charge. R36 is near to the α-phosphate of AD(M)P along the entire catalytic reaction, stabilizing both the transition state and products. Overall, interactions between these arginine side chains and several backbone amides and the Mg²⁺ with the phosphates of the substrate make up a well-organized, asymmetric active site enabling efficient reversible P-transfer with a wide TSE (Figure 3g).

Umbrella sampling and committor analysis buttress wide TSE for Adk with Mg²⁺

To verify our computational finding of a wide TSE in the fully assembled enzyme with Mg²⁺ in contrast to a narrower TSE without a divalent metal, we performed QM/MM umbrella sampling and a committor analysis as described in methods. The umbrella sampling free energy profiles (Figure 4a-d) corroborate the key difference in the range of the TSE. Noteworthy, the activation free energies are smaller for the umbrella sampling when compared to the MSMD. The largest difference being that for the reaction in the absence of Mg(Figure 1c). Higher barriers when comparing the results using Jarzynski’s to those obtained with umbrella sampling is not totally unexpected. It is well known that exponential averaging tends to overestimate free energy barriers, and particularly when they are too high (Park et al., 2003).

Figure 4:

Notably, the calculated difference in activation barriers with and without Mg from these umbrella simulations of about 9.5 kcal/mol are in good agreement with the experimentally determined ones of ≥ 11 kcal/mol, (P-transfer with Mg is > 500 s⁻¹ (Kerns et al., 2015) compared to 7.5×10⁻⁴ s⁻¹ measured here).

A consequent comparative commitment analysis that calculates the probability of reaching either reactants or products displays a very shallow change along the reaction coordinate in the presence of Mg²⁺, further corroborating a wide TSE for the P-transfer step in the fully active enzyme, in stark contrast to a narrow and steep change in the committors indicative of a narrower TSE without Mg²⁺ (Figure 4d). The TSE is defined in function of the committor distribution (Figure 3 — supplement 9), showing the TSE is at −0.1 Å of the reaction coordinate for the reaction without Mg²⁺. In contrast, the TSE is wider for the reaction with Mg²⁺ (reaction coordinate values in the range from 0.0 to 0.2 Å (Figure 3 — supplement 9 and Table 1).

Experimental characterization of the activation parameters of the P-transfer step

We felt the need to additionally experimentally test our major new finding from the QM/MM simulations: a delocalized TSE for the fully active enzyme. This feature would result in a lowering of the entropic barrier for the chemical step, thereby contributing to the enzyme-catalyzed rate enhancement. To avoid the issues of accuracy of QM/MM simulations that are well documented (Acevedo & Jorgensen, 2010; Elstner, 2007; Gaus et al., 2014; Roston et al., 2018), our system has the advantage of comparing the exact same reaction coordinate with the only difference being the presence/absence of a single atom, Mg²⁺, thereby studying differences rather than absolute values. Since the reaction with Mg²⁺ (fully active enzyme) revealed a more pronounced delocalization than the reaction in its absence, the divalent cation is predicted to lower the activation entropy. To test this prediction, we experimentally determined the enthalpic and entropic contributions to the chemical reaction barrier in Adk by measuring its temperature dependence in the presence and absence of Ca²⁺. The experimental trick of replacing Mg²⁺ with Ca²⁺ was used since the chemical step with Mg²⁺ is too fast to be experimentally measured, and previous studies showed that Ca²⁺ is an appropriate mimic to selectively probe the chemical step (Kerns et al., 2015).

The resulting Eyring plots for the two experiments (Figure 5a) deliver the enthalpic (ΔH^‡ of 15.3 ± 0.5 and 16.7 ± 1 kcal/mol) and entropic (ΔS^‡ = −-5.7 ± 1.6 and −19.1 ± 3.1 cal/mol/K) contributions to the chemical reaction barrier in the presence and absence of Ca²⁺, respectively. The presence of Ca²⁺ reduces both barriers considerably. The decrease in the enthalpic barrier by Ca²⁺, which is even more pronounced with the optimal Mg²⁺ ion, is not new and has been well-established in the literature for P-transfer reactions, including the Adk chemical reaction_(Kerns et al., 2015; Pérez-Gallegos et al., 2017; Rosta et al., 2014; L. Yang et al., 2012). The striking new result, central to this work, is the large decrease in the entropic barrier in the presence of a divalent cation, thereby strongly supporting our key finding of a broader TSE from the QM/MM simulations.

Figure 5:

We were able to further test our QM/MM simulations with a second designed set of experiments by measuring the pH dependence of the P-transfer step. We find that the reaction rate increases with higher pH in the presence of Ca²⁺, whereas without a metal the opposite trend is observed (Figure 5b). These experimental results match the simulation results: the fully-charged nucleotides state was the most reactive in the presence of the metal, whereas the monoprotonated state had a much higher free-energy barrier. In contrast, without metal only the monoprotonated nucleotides were reactive. Notably, at low pH (pH 5), the P-transfer step becomes rate limiting even with Mg²⁺, allowing for an estimate of the rate enhancement of Mg²⁺ versus Ca²⁺ of 24-fold.

Finally, we measured the forward and backward reaction rate constants of several arginine mutants that were identified as important for the P-transfer from the QM/MM simulations: R36K, R85K, R124K, R150K, and R161K (Table 2). The large decrease in the catalytic rate for each single Arg to Lys mutation further highlights the notion of a highly-choreographed active site for a well-coordinated P-transfer (Figure 3g), in which arginine residues that are not even directly coordinated to the transferring phosphate play an equally important role.

Table 2.

Discussion

Adk is an impressively efficient enzyme, being able to accelerate the P-transfer reaction over a billion times. Previous work showed the complexity of the Adk free-energy landscape for this two-substrate, two-product reaction (Figure 1a), providing direct experimental evidence for the widely accepted concept of a multidimension free energy landscape with multiple intermediates and transitions states for the full enzymatic cycle, already described and demonstrated by Benkovic and Hammes-Schiffer (Benkovic et al., 2008). Full lid-closure results in an excellent preorientation of the donor and acceptor groups of the two bound nucleotides and the aiding active site residues in the enzyme/substrate complex, in agreement with the general view of necessary pre-organization in enzyme catalysis. Lid-opening had been identified as the rate-limiting step in the enzymatic cycle (Kerns et al., 2015), and the corresponding TSE for this conformational change (Stiller et al., 2019) and the high-resolution structure of next intermediate, the EP-complex before product dissociation, recently been described (Stiller et al., 2022).

For a comprehensive understanding of the complete free-energy landscape of Adk catalysis, we here investigate the key step, catalyzing the actual chemical step of P-transfer. While ES and EP complexes such as Adk bound to ADP/Mg²⁺ can be structurally characterized by traditional experimental methods as they represent a minimum in the free-energy landscape, the reaction path from substrate to product involving breaking and forming covalent bonds (chemical step) and traversing the crucial transition state can only be “visualized” by Quantum mechanics based molecular simulations. The power of such simulations to examine the chemical steps in enzymes, including P-transfer reactions, has been extensively documented (Jin et al., 2017; Lai & Cui, 2020b; Mokrushina et al., 2020; Pérez-Gallegos et al., 2015, 2017; Roston & Cui, 2016; Valiev et al., 2007). The focus in the literature has been on the comparison of the enzyme-catalyzed and uncatalyzed reaction with respect to the bond character at the transitions state (i.e., associate versus dissociative or tight versus loose) (Admiraal & Herschlag, 1995; Hengge, 2002; Kamerlin et al., 2013; Lassila et al., 2011). QM calculations for uncatalyzed reactions have been well documented in the literature showing narrow and symmetric TSE (Klähn et al 2006: Wang et al. 2015).

Here, we uncover a central new result that provokes a modified transition-state theory. Enzymes, due to their macromolecular nature, provide a fundamentally different, advantageous way to catalyze these chemical reactions compared to the uncatalyzed reaction by employing a broad TSE. Many different molecular configurations can be accommodated in the transition-state region with comparable energies via collective motions, spanning a 1 Å range along the reaction coordinate for the transferring phosphate that travels a total distance of only 2 Å during the P-transfer. This features resemblance the now well-established conformational sampling of proteins in ground states such as ES complexes. Our findings explain why enzymes do not face an entropic bottleneck for catalysis. Furthermore, as enzyme active sites are asymmetric in contrast to the symmetric nature of the solvent for uncatalyzed reactions, we find that the TSE comprises also highly asymmetric conformations. Our findings help to resolve the controversy about the nature of the transition state in enzyme-catalysed P-transfer reactions between theory and experiments (Lassila et al., 2011). The complex nature of the active site of enzymes, in contrast to simple solvent, results in different mechanisms in the enzyme catalyzed reaction.

We note that a previous QM/MM minimum energy path calculations for a different Adk enzyme, E.coli Adk, using a semiempirical method had proposed a different mechanism with a stable metaphosphate intermediate, reporting an even lower energy for this metaphosphate intermediate than the ES- and EP-complexes (Shibanuma et al., 2020). In complex systems, such as enzymes, it is possible to observe artificial local minima when using minimum energy path searching strategies due to inadequate sampling (Mendieta-Moreno et al., 2015; Quesne et al., 2016). From our experimental NMR and x-ray data we know that such a stable metaphosphate does not exist in Adk-catalyzed reactions, highlighting the importance of experimental verification of simulations as performed here, and the use of extensive sampling with proper thermodynamic treatment

More recently, another group performed QM/MM simulations for E.coli Adk using the semiempirical AM1/d-phoT method (Nam et al. 2007) and umbrella sampling, and measured and computed rate reduction in the chemical step by similar Arg mutants as reported here (Ojeda-May et al., 2021 and Dulko-Smith et al., 2023). In contrast to the paper described above (Shibanuma et al., 2020), their results are in agreement with our data. First, the experimentally measured rate reductions in (Ojeda-May et al., 2021 and Dulko-Smith et al., 2023) are in full agreement with our results. Second, our activation barriers from umbrella sampling are qualitatively similar to their values, with the obvious differences in the enzyme species (we study the thermophilic Aquifex Adk), and their choice of analyzing the QM/MM simulations in the opposite direction ATP/AMP to 2 ADP molecules. The latter lead them to singly protonate one oxygen on ATP, since this “back reaction” is faster for the mono-protonated state, in agreement with our simulations. This nice agreement in both experiments and computation by our two groups strengthens the validity of reported results.

Intriguingly, the Nam/Wolf-Watz team (Ojeda-May et al., 2021 and Dulko-Smith et al., 2023) and our group focus on completely different interpretations of the data! Nam/Wolf-Watz propose an interconnection, synchronization, between the chemical step and lid-opening. This language can be misleading since it can imply that the chemical step directly effects the opening, which of course cannot be the case. There is no “memory” of the chemical step, the energy from the closing/opening step is dissipated within picoseconds. Lid-opening/closing is an energetically fully independent step from the P-transfer step, see our Figure 1a. We agree with the obvious result/interpretation by the authors that the Arg residues in the active site are important for both the P-transfer and the consequent lid opening, but there is no synchronization between these independent steps. In contrast, we discover here a fundamental concept for rate enhancement by an optimal enzyme, the reduction in the activation entropy by a wide TSE. New experiments were triggered by our finding that then delivered experimental validation of this concept.

A wide TSE agrees with the definition of the TS being a surface on the potential energy surface instead of a saddle as elegantly discussed in a theoretical paper by Nussinov and colleagues (Ma et al., 2000). Our results deliver concrete evidence for this fundamental concept from QM/MM calculations and commitment analysis on Adk that were further buttressed by the consequently designed experiments. Although we have only quantitatively demonstrated this TSE concept here for Adk catalyzed P-transfer, we hypothesize such mechanism to be more general for enzyme catalysis as it is rooted in the macroscopic nature of proteins. A conformationally spread-out TSE is likely one reason why enzymes are much bigger than a few active-site residues: On one hand residues remote from the site of the chemical action allow for conformational flexibility without the protein unfolding, on the other hand they allow for increasing the probability of highly-choreographed motions along the reaction coordinate (Hammes-Schiffer & Benkovic, 2006; Klinman & Kohen, 2013; Saen-Oon et al., 2008; Schramm & Schwartz, 2018) while reducing excursions into unproductive conformations (Otten et al., 2020). The latter mechanism is directly demonstrated here in the TSE, as the configurations seen in the simulations are along the P-transfer reaction coordinate. We note that we focus here on equilibrium effects and an atomistic description of the TSE, this is not related to the fundamental work on dynamical effects manifested in transmission coefficients due to barrier recrossing (Antoniou & Schwartz, 2016; Zinovjev & Tuñón, 2017). The latter has been quantitatively demonstrated to play a minor role.

A question often being asked is why would nature evolve a conformational step to be rate limiting, as the lid-opening step for Adk? We propose that this is the “price” the enzyme pays for stabilizing the ES complex to enable efficient rate acceleration of the chemical step, and for suppressing detrimental alternate reactions such as phosphate hydrolysis. Release of product requires disassembly of such closed, low-energy structure. These principles have now been directly revealed on an atomic scale for Adk: extensive electrostatic interactions of five Arg and one Lys with the phosphoryl groups to efficiently accelerate the P-transfer step, which then need to be broken for lid opening to allow product release. The effect of different Arg mutations on the lid opening has been recently nicely demonstrated (Ojeda-May et al., 2021 and Dulko-Smith et al., 2023).

As long as the rate of this conformational transition (k_open) and hence k_cat is not limiting for organismal fitness, there is no evolutionary pressure to further enhance it.

The importance of TSE including multiple transition pathways for protein folding (Frauenfelder et al., 2006; Royer, 2008), and more recently for protein conformational transitions within the folded state (Pontiggia et al., 2015; Stiller et al., 2019; Tsai & Nussinov, 2014) has been well documented. TSE in enzyme catalysis of chemical steps where covalent bonds are broken and formed as demonstrated here further promotes a unifying concept of protein folding and function, as the same principal concept of TSE is seen for these different types of conformational changes. Consequently, protein folding and function is embedded in a unifying energy landscape which is simply altered by ligands and solvent.

Materials and Methods

Computational Set up of the system

The crystal structure of adenylate kinase from Aquifex aeolicus in complex with Ap5A coordinated to Zn²⁺ (PDB ID 2RGX (Henzler-Wildman et al., 2007)) was used as model. The Zn²⁺ ion was replaced by Mg²⁺ and the ADPs structures were built by modifying the structure of Ap5A, keeping the crystallographic coordination of the cation (2 oxygens from ADPs and 4 waters). For the simulation in low pH and without the magnesium ion, the nucleophilic oxygen of AD(T)P was protonated. Standard protonation states were assigned to all titratable residues (aspartate and glutamate residues are negatively charged, lysines and arginines residues are positively charged). Histidine protonation was assigned favoring formation of hydrogen bonds in the crystal structure. Each protein was immersed in a truncated octahedral box of TIP3P water (Jorgensen et al., 1983) and neutralized with Na⁺ ions (3 to 5 ions, depending on the system).

Classical Simulation Parameters

The parameters used for all protein residues were those corresponding to AMBER force field ff99sb (Hornak et al., 2006). For the non-protein residues, such as ADP and Mg²⁺, parameters were taken from refs. Meagher et al. (Meagher et al., 2003) and Allnér et al. (Allnér et al., 2012). All simulations were performed with periodic boundary conditions. Pressure and temperature were kept constant with the Berendsen barostat and thermostat, respectively (Berendsen et al., 1984). The SHAKE algorithm (Ryckaert et al. 1977) was used to keep Hydrogen atoms at their equilibrium bond distances, and a 2 fs time step was employed, except for the system without Mg²⁺ for which it was set to 1 fs. The ADP-ADP system was initially subjected to a classical equilibration protocol, followed by QM/MM equilibration and production simulations. Starting from the classically equilibrated ADP-ADP, the mono-protonated systems with ADPH-ADP was built, by adding the H+ and removing a Na⁺ from the solvent, and subsequently the system without Mg²⁺ was built by removing the ion and adding two Na⁺ to the solvent. Each of these systems was again equilibrated classically and then used for forward QM/MM simulations. The backward simulations were equilibrated just using QM/MM simulations because the starting point was the product of the forward SMD simulation with the lowest work at product RC. Classical equilibration followed consisted of: first, 2000 steps of conjugate gradient optimization, followed by 100ps MD simulation at NVT condition, where the system’s temperature was slowly raised to 300 K; The third step consisted of 100ps MD simulation at NPT conditions to equilibrate the system’s density. During the temperature and density equilibration process, the protein α-carbon atoms were restrained by a 1 kcal/mol harmonic potential. Finally, non-restrained classical simulation was performed at least for 1 ns at NPT conditions to obtain an equilibrated structure to be used as the QM/MM input.

QM/MM simulations

were performed using a 1 fs timestep without SHAKE (Ryckaert et al. 1977). To equilibrate the system with the QM/MM hamiltonian, we performed first a conjugate gradient QM/MM optimization, followed by a 50 ps QM/MM MD equilibration process at NVT condition, subsequently each production simulation was performed at NVT condition (at least 100 ps saving coordinates every 1 ps). QM/MM calculations were performed using the additive scheme implemented in the sander module of the AMBER program suite (Case et al., 2005). The quantum (QM) region was treated with the self-consistent charge-density functional tight-binding (SCC-DFTB) approach (Elstner, 2007), usually considered a semi-empirical QM method, as implemented in Sander (Seabra et al., 2007; Walker et al., 2008). The SCC-DFTB parameter set involving carbon, hydrogen, oxygen, nitrogen, phosphate, and magnesium was employed throughout this study (Cai et al., 2007; Y. Yang et al., 2008). All the other atoms were described with AMBER ff99SB force field (Hornak et al., 2006). The interface between the QM and MM portions of the system was treated with the scaled position-link atom method. The electrostatic interactions between the QM and MM regions were treated using electrostatic embedding using a simple mulliken charge - RESP charge interaction. In the case of the van der Waals interactions, these were calculated classically.

The QM region of the complete system consists of the diphosphate moiety of both ADP molecules, one Mg2+ ion, and the four water molecules that coordinate the metal ion in the crystal structure (Figure 1 – supplement 1).

Free energy determination strategy using MSMD and Jarzynski’s equation

The free- energy profiles were constructed by performing constant velocity multiple steered molecular dynamics (MSMD) simulations (i.e., stiff string approximation), and using Jarzynski’s equality (Jarzynski, 1997), which relates equilibrium free-energy values with the irreversible work performed over the system along a user defined reaction coordinate that drives the system from reactants to products. In the present study, the reaction coordinate (ξ) was chosen as:

Calculations were performed using a force constant of 300 kcal/mol/Å and pulling velocities of 0.05 Å ps⁻¹. To reconstruct the free-energy profile of the phosphate transfer reaction, two sets of at least 10 SMD runs were performed starting from equilibrated QM/MM MD structures corresponding to either the system in 1) the Adk/ADP-ADP state (presently defined as forward reaction), and 2) the Adk/ATP-AMP state (backward reaction). Each set of work profiles was used to obtain the corresponding forward and back free-energy profiles with Jarzynski’s equality (Jarzynski, 1997), and finally the two profiles were combined to obtain the complete curve. Since when using MSMD in combination with Jarzynski’s inequality, the system starts from equilibrium conditions and is driven along the RC further from equilibrium as the reaction proceeds, and therefore the FEP estimate increasingly overestimates the real FEP. The best way to combine forward and back reaction is performed by keeping the initial segment (and thus lower) FEP of each forward and back reaction estimates. This strategy was successfully used to obtain enzymatic reaction free-energy profiles in previous works (Crespo et al., 2005; Defelipe et al., 2015).

Along the present work, the ADP in the ATP-lid will be called as AD(T)P and the other, in the AMP-lid, as AD(M)P. The reaction involves the two P-O bonds, one forming and the another breaking, the two bonds sharing the same phosphorus atom from the transferring phosphate (Figure 1b). In the former, an oxygen atom of the AD(T)P reacts and we will name it as O_attacking. The breaking bond is described by the oxygen atom for the AD(M)P and it will be called O_leaving.

Higher level DFT(PBE) free energy calculations

The free-energy profile of the ADK with Mg²⁺ was also computed using a higher level of theory. Here, the QM system was described with PBE functional (Perdew et al., 1996) using a DZVP basis set as implemented in the GPU based code LIO that works with amber (Nitsche et al., 2014). To take advantage of previous DFTB sampling, five lowest work vs RC profiles of the forward and back reaction were selected, and segmented in 10 windows or stages. For each window, the free energy profile at the DFT(PBE) level was computed using SMD with a steering velocity of 0.05 Å/ps. Windows were combined to obtain both forward and back FEP and the final profile was obtained by joining them. Details of the staged strategy can be found in references Ozer et al. (Ozer et al., 2010, 2012) and Ganguly et al. (Ganguly et al., 2020).

Commitment analysis

To verify the true nature of the observed TS we performed a commitment analysis. We first selected 20 structures for each of 11 points along the RC in the TS region (every 0.1 Å from −0.6 to 0.6 Å on the RC). Structures were equilibrated at the corresponding RC position for 5 ps using a harmonic restrain potential of 300 kcal/mol/Å. From the final equilibrated structure we used the modified Berendsen thermostat (Bussi et al, 2007) to assign random initial velocities, and performed 10 non restrained simulations. Each simulation was followed and determined whether it reached the reactants or products. Finally, the resulting probability (estimated from the observed frequency) to reach product for each point along the RC was computed.

Umbrella Sampling simulations

Finally, to avoid possible bias due to the use of MSMD to obtain the FEP, we also computed the corresponding profile using Umbrella Sampling (US). The initial structures for each window were randomly selected from the lowest work MSMD trajectories along the RC. The force constant for the harmonic restraint was set to 300 kcal/mol/Ų, ensuring sufficient confinement without inducing artificial behavior. Windows were separated every 0.1 Šalong the whole RC. Each window simulation started with a equilibration run of 60 ps at NPT condition, followed by sampling over a production run of 40 ps, with coordinates recorded at 1 ps intervals and the reaction coordinate values every 1 fs. Simulations were performed at the previously described SCC-DFTB level of theory. The weighted histogram analysis method (WHAM) was employed to reconstruct the free energy landscape from the individual window simulations (Kumar et al. 1992; Grossfield, A.).

Principal Component Analysis (PCA)

One single trajectory containing the TSE structures of MSMD and Umbrella Sampling was created by aligning the atoms of the P-O_leaving and P-O_attacking to a reference. The reference was the reactant structure (ADP/ADP). This was done in order to compare the TSEs using the same coordinate system. The PCA was calculated over the transferring phosphate atoms, Mg2+ ion, P_α and O3α atoms from AD(M)P and P_β and O3β atoms from AD(T)P. Those atoms were selected because they are directly involved in the chemical reaction.

Steady-state kinetics measurements

Steady-state kinetics measurements for Aquifex aeolicus Adk (Adk) were performed at different pH values in the absence of divalent metals (using 50 mM of EDTA) and the presence of Mg²⁺ or Ca²⁺. In all cases the reaction was started by the addition of 8 mM ADP and saturating metal concentrations were used if applicable (32 mM at pH 5 and 8 mM at pH 7 and 9, respectively). The enzyme concentration varied between 30 nM and 200 μM depending on the rate of interconversion; measurements were collected at 25 °C. The amount of product was quantified with high-pressure liquid chromatography (HPLC). Protein precipitated by quench (30% TCA + 6 M HCl mixture) was separated with Spin-X centrifugal tube filters (Costar), filtered supernatant was diluted to avoid HPLC detector saturation, and the pH was brought to 6.0 to achieve optimal separation. The samples were analyzed on an HPLC system (Agilent Infinity 1260) with a high-precision autosampler (injection error <0.1 μL) and analytical HPLC column ACE (i.d. 2.4 mm, length 250 mm, C18-AR, 5 Å pore size) and separated with isocratic elution with potassium phosphate mobile phase (100 mM, pH 6.1). The observed rate constants were determined from 8 – 15 data points for each temperature and/or pH value using initial rate analysis. The values and uncertainties (s.d.) shown in Figure 5 and Table 2 were determined from least-squares linear regression. Similarly, values and uncertainties (s.d.) of ΔH^‡ and ΔS^‡ were extracted from linear regression of the data points presented in Figure 5a.

Temperature Dependence of catalysis and Mutant activity

Experiments were essentially performed as described previously⁶ and above. In short, the steady-state kinetics measurements were collected at temperatures between 20 and 60 °C for the temperature-dependency and at 25 °C for the mutants. For the temperature dependence, the samples contained 4 mM ADP and equimolar (with nucleotide) concentrations of calcium or 50 mM EDTA. The enzyme concentration was varied between 50 nM – 1 μM (with Ca²⁺) or 100 – 200 μM (with EDTA); buffer was 100 mM HEPES, pH 7.0, and 80 mM KCl. The amount of product produced over time (10-16 minutes) was quantified with HPLC and observed rates were extracted as described above. A similar approach was used for the Arg-to-Lys Aadk mutants, where the reaction was measured in both directions starting with either 4 mM ADP or ATP/AMP and equimolar (with nucleotide) concentrations of Mg²⁺.

Supplementary Figures

Figure 1 – supplement 1.

Figure 1 – supplement 2.

Figure 1 – supplement 3.

Figure 1 – supplement 4.

Figure 1 – supplement 5.

Figure 1 – supplement 6.

Figure 3 – supplement 1.

Figure 3 – supplement 2.

Figure 3 – supplement 3.

Figure 3 – supplement 4.

Figure 3 – supplement 5.

Figure 3 – supplement 6.

Figure 3 – supplement 7.

Figure 3 – supplement 8.

Figure 3 – supplement 9.

Acknowledgements

This work was supported by the Howard Hughes Medical Institute (HHMI) to D.K. We thank to “High-Performance Computing Center” (CeCAR, https://cecar.fcen.uba.ar/) of the Facultad de Ciencias Exactas y Naturales at the University of Buenos Aires for the computational resources. We are grateful to the Brazilian Biosciences National Laboratory (LNBio), part of the Brazilian Center for Research in Energy and Materials (CNPEM) for accessibility to the High-Performance Computing Cluster and scientific infrastructure.