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). Ironically, the piece missing to understand the enzyme’s catalytic power is the mechanism of the actual P-transfer step, which would take about 7000 years without the enzyme (Stockbridge & Wolfenden, 2009) compared to > than 5,000 s^-1 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). Ten 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²⁺.

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 and 1e). 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 2:

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 1 – supplements 6-8). 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.

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 3-5). 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).

Stimulated by experts in the field that keep pushing the methodology, 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 9). The results are in excellent agreement with the free energy profiles obtained before (Figures 1a,e and Figures 2a,c) thus further supporting the computational results.

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

We felt the need to experimentally test our major new finding from the QM/MM simulations: a delocalized TSE. 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 4a) deliver the enthalpic (ΔH^# of 64.1 ± 2.2 and 70.0 ± 4.1 kJ/mol) and entropic (ΔS^# = -23.8 ± 6.8 and -80.0 ± 12.8 J/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 4:

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. The reaction rate increases with higher pH in the presence of Ca²⁺, whereas without a metal the opposite trend is observed (Figure 4b). 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 land-scape 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 well 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 next intermediate, the EP complex before product dissociation (Stiller et al., 2022), recently been described.

However, for a comprehensive understanding of the complete free-energy landscape of Adk catalysis, the key step, catalyzing the actual chemical step of P-transfer, had been missing. 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). 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 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 further note that a previous QM/MM minimum energy path calculations for Adk using a semiempirical method had proposed a different mechanism with a stable meta-phosphate intermediate with 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 NMR and x-ray experimental data we know that such stable meta-phosphate 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.

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 the free-energy profiles of the QM/MM calculations 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 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 5 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. However, 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 now 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. Each protein is uniquely described by the corresponding free energy landscape that includes folded and unfolded conformational substates and is defined by the amino acid sequence and the surrounding solvent. 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 Aq-uifex 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 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). The parameters used for the protein were those corresponding to AMBER force field ff99sb (Hornak et al., 2006). For the non-protein residues, as ADP and Mg²⁺, the parameters were taken from refs. Meagher et al. (Meagher et al., 2003) and Allnér et al. (Allnér et al., 2012).

Simulation Parameters

(a) Classical force field and Scheme: All simulations were performed with periodic boundary conditions; the SHAKE algorithm was used to retain the hydrogen atoms at equilibrium bond lengths, allowing the use of a 2 fs time step. Pressure and temperature were kept constant with the Berendsen barostat and thermostat, respectively (Berendsen et al., 1984). All the simulations were equilibrated following the same protocol: 1) 2000 steps of optimization using a conjugate gradient algorithm; 2) 100 ps MD simulation at NVT condition where the system’s temperature was slowly raised to 300 K; 3) 100 ps MD simulation at NPT conditions to equilibrate the system’s density. The protein α-carbon atoms were restrained by a 10 kcal mol^-1 harmonic potential for the thermalization and by a 1 kcal mol^-1 harmonic potential for the density equilibration. Afterwards, a production simulation was performed for 7.5 ns at NPT conditions. (b) QM/MM parameters: QM/MM simulations were performed starting from selected snaphots from the classical production simulation. 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, subsequently each production simulation was performed. 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 QM region of the complete system consisted of the diphosphate moiety of both ADP molecules, one Mg²⁺ ion, and the four water molecules that coordinate the metal ion in the crystal structure (Figure 1 – supplement 1).

Free energy determination strategy

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^-1 Å^-1 and pulling velocities of 0.05 Å ps^-1. To reconstruct the free-energy profile of the phosphate transfer reaction, two sets of 10 SMD runs were performed starting from equilibrated QM/MM MD structures corresponding to either the system in 1) the Adk/ADP-ADP state (for forward reaction), and 2) the Adk/ATP-AMP state (for 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. This strategy was successfully used to obtain enzymatic reaction free-energy profiles in previous works (Crespo et al., 2005; Defelipe et al., 2015). 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).

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 4 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 4a.

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²⁺.

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.

Competing interests

D.K. is co-founder of Relay Therapeutics and MOMA Therapeutics. The remaining authors declare no competing interests.

Supplementary Figures

Figure 1 – supplement 1.

Figure 1 – supplement 2.

Figure 1 – supplement 3.

Figure 1 – supplement 4.

Figure 1 – supplement 5.

Figure 3 – supplement 1.

Figure 3 – supplement 2.

Figure 1 – supplement 6.

Figure 1 – supplement 7.

Figure 1 – supplement 8.

Figure 3 – supplement 3.

Figure 3 – supplement 4.

Figure 3 – supplement 5.

Figure 1 – supplement 9.