The role of the Met20 loop in the hydride transfer in Escherichia coli dihydrofolate reductase

A key question concerning the catalytic cycle of Escherichia coli dihydrofolate reductase (ecDHFR) is whether the Met20 loop is dynamically coupled to the chemical step during catalysis. A more basic, yet unanswered question is whether the Met20 loop adopts a closed conformation during the chemical hydride transfer step. To examine the most likely conformation of the Met20 loop during the chemical step, we studied the hydride transfer in wild type (WT) ecDHFR using hybrid quantum mechanics-molecular mechanics free energy simulations with the Met20 loop in a closed and disordered conformation. Additionally, we investigated three mutant forms (I14X; X = Val, Ala, Gly) of the enzyme that have increased active site flexibility and donor-acceptor distance dynamics in closed and disordered Met20 loop states. We found that the conformation of the Met20 loop has a dramatic effect on the ordering of active site hydration, although the Met20 loop conformation only has a moderate effect on the hydride transfer rate and donor-acceptor distance dynamics. Finally, we evaluated the pKa of the substrate N5 position in closed and disordered Met20 loop states and found a strong correlation between N5 basicity and the conformation of the Met20 loop.

A key question concerning the catalytic cycle of Escherichia coli dihydrofolate reductase (ecDHFR) is whether the Met 20 loop is dynamically coupled to the chemical step during catalysis. A more basic, yet unanswered question is whether the Met 20 loop adopts a closed conformation during the chemical hydride transfer step. To examine the most likely conformation of the Met 20 loop during the chemical step, we studied the hydride transfer in wild type (WT) ecDHFR using hybrid quantum mechanics-molecular mechanics free energy simulations with the Met 20 loop in a closed and disordered conformation. Additionally, we investigated three mutant forms (I14X; X ‫؍‬ Val, Ala, Gly) of the enzyme that have increased active site flexibility and donor-acceptor distance dynamics in closed and disordered Met 20 loop states. We found that the conformation of the Met 20 loop has a dramatic effect on the ordering of active site hydration, although the Met 20 loop conformation only has a moderate effect on the hydride transfer rate and donor-acceptor distance dynamics. Finally, we evaluated the pK a of the substrate N5 position in closed and disordered Met 20 loop states and found a strong correlation between N5 basicity and the conformation of the Met 20 loop.
Escherichia coli dihydrofolate reductase (ecDHFR) 3 serves as a benchmark enzyme for the study of static and dynamic effects in biocatalysis (1). DHFR catalyzes the reduction of dihydrofolate (DHF or H 2 folate) by nicotinamide adenine dinucleotide phosphate hydride (NADPH) to form tetrahydrofolate (H 4 folate) and NADP ϩ (Scheme 1). Its principal function is to maintain intracellular pools of H 4 folate, a compound essential for the biosynthesis of purines, thymine nucleotides, and several amino acids.
The Met 20 loop in ecDHFR is a flexible structural motif that closes and opens the active site cleft as ligands bind or dissociate from the enzyme. The Met 20 loop dynamics have been suggested to occur on the pico-to millisecond timescale (2), and an important question is whether any of its motions are coupled to catalysis (3)(4)(5). It has been proposed that the Met 20 loop (residues 9 -23) may adopt three principle states, open, closed, and occluded (6), and the conformation of the Met 20 loop has been proposed to correlate with catalysis (6 -11). For instance, Brooks and co-workers (12)(13)(14) suggested that rapid hydride transfer is correlated with a selected subset of Met 20 conformational states. This raises the question of whether the Met 20 loop is closed during the chemical step and more generally whether the Met 20 loop plays a role during hydride transfer.
To investigate this question, we studied the hydride transfer in wild type (WT) ecDHFR and three Met 20 loop mutant enzymes, namely I14V, I14A, and I14G. These mutants have been shown to possess increased active site flexibility relative to the WT enzyme and, in particular, enhanced donor-acceptor distance (DAD) dynamics (15,16). This DAD dynamics perturbation has a rather radical effect on the temperature dependence of the kinetic isotope effect (KIE) (15,17). In the WT enzyme, the KIE is temperature-independent, whereas the KIE becomes increasingly temperature-dependent as the steric effect of the mutation becomes more pronounced. In the crystal structures of these mutant enzymes, the Met 20 loop adopts an increasingly open and disordered conformation as a function of the mutation (Fig. 1) (16). The analogous ternary complexes for all three mutations of Ile 14 in ecDHFR were crystallized, and their structures were determined (16). In the I14A and I14V mutants, the structures were obtained in the P2 1 space group where the Met 20 loop is observed in a disordered conformation as this conformation is seemingly stabilized by crystal contacts. These structures are rather similar to the analogous WT structures of Sawaya and Kraut (6) in the C2 and P2 1 space groups. In contrast, the crystal structure of I14G ecDHFR was crystallized in the P2 1 2 1 2 space group where both the Met 20 loop and the F-G loop, which is in direct contact with the Met 20 loop (Fig. 1), are partially disordered. It is possible that a transient population of the open conformational state plays a role in the catalytic cycle because the open conformation enables both substrate binding and product release (3, 18 -20). However, replacement of a part of the Met 20 loop (residues 16 -19) with a single glycine residue still yields an enzyme that exhibits some catalytic activity (7), suggesting that the loop is not essential for basic catalytic activity. In our previous work, the above mentioned mutant crystal structures with a disordered Met 20 loop were not used for modeling purposes as we focused on the closed conformation. Indeed, it is not clear whether the Met 20 loop plays a role during the hydride transfer step in general (7) and whether it influences the DAD in particular.
In this work, we present free energy profiles along with structural and dynamic analyses for the hydride transfer in WT ecDHFR and three of its mutants (I14V, I14A, and I14G). The WT enzyme is treated in a closed conformation as is observed in the ternary complex crystal structure as well as with a disordered Met 20 loop conformation (6). The mutant enzymes are also constructed in both the closed and disordered Met 20 loop conformations. To quantify the effect of conformation on chemistry, we used a hybrid quantum mechanics (QM)/molec-ular mechanics (MM) Hamiltonian especially designed for the DHFR reaction (21). The free energy profiles were obtained from molecular dynamics (MD) umbrella sampling simulations and were compared with available experimental pre-steadystate single-turnover rates (15). Finally, we estimated the pK a for the N5 position of the folate moiety with the Met 20 loop in the closed and disordered conformations.

Free energy profiles
We obtained the classical mechanics potential of mean force (PMF) profiles at 25°C for the catalyzed hydride transfer from NADPH to H 3 folate ϩ (N5-protonated DHF) in the WT ecDHFR and three I14X mutant ecDHFRs, namely I14V, I14A, and I14G. In Fig. 2A, we present the PMFs for the WT enzyme as well as for the I14X mutant enzymes, all with the Met 20 loop in the closed conformation. The classical free energy of activation for the WT enzyme, ⌬G ‡ , is 17.4 Ϯ 0.3 kcal/mol, whereas the reaction free energy, ⌬G r , is Ϫ4.5 Ϯ 0.3 kcal/mol. The barrier obtained is similar to that obtained in our previous study (16.8 kcal/mol), although the reaction is less exergonic than previously observed (Ϫ9.7 kcal/mol) (16,21). Possible sources for the slight differences are the different initial crystal structures used in the two studies, higher ionic concentration in the

Role of the Met 20 loop in the hydride transfer in E. coli DHFR
current work, and statistical uncertainty. We stress that addition of nuclear quantum effects on the hydride transfer typically reduces the barrier for this reaction by 2-3 kcal/mol (21,22), bringing the computed value into agreement with the phenomenological value estimated from single-turnover rates and the Eyring equation (13.4 -14.3 kcal/mol) (15,(23)(24)(25). We have previously shown that the nuclear quantum effects are similar for WT and the I14X mutants (16,22). The presence of a slightly less bulky side chain in the I14V ecDHFR system causes an increase in the free energy barrier (18.3 Ϯ 0.4 kcal/mol) relative to WT ecDHFR in mutant simulations with the Met 20 loop closed (Figs. 1A and 2A). Further truncation of the side chain brings about even higher free energy barriers as evident by the ⌬G ‡ value for the I14A and I14G systems (18.4 Ϯ 0.4 and 18.8 Ϯ 0.4 kcal/mol, respectively). The effect of the mutations is similar to our earlier work and is generally in accord with the experimental trends, although the effect of mutation to Val is slightly overestimated ( Table 1) (15)(16)(17). We note that the experimental single-turnover rates are pH-dependent and include conformationally induced fit during DHF binding, the protonation at N5, and the hydride transfer step, whereas the current simulations only account for the hydride transfer step. In Fig. 2B, we show an analogous plot but with the Met 20 loop in a disordered conformation (Fig. 1B). A similar trend in the PMFs is seen in these profiles (18.3 Ϯ 0.4, 18.5 Ϯ 0.1, 18.8 Ϯ 0.4, and 20.9 Ϯ 0.4 kcal/mol for WT*, I14V*, I14A*, and I14G*, respectively). The free energy barriers are somewhat higher for all models with the Met 20 loop in a disordered conformation than with the Met 20 loop in the ordered conformation (see supplemental Fig. S3). We ascribe this to the greater reorganization cost in climbing from the ground state (GS) to the transition state (TS) required for the more disordered structures. We note that the limited accuracy in our computed free energy barriers (approximately Ϯ0.4 kcal/mol) precludes us from drawing quantitative conclusions regarding possible correlation between the effects of Met 20 loop conformation and I14X mutations, such as additivity or cooperativity.

Structural analysis
Conformational substates-We examined different loop conformations that are sampled in each of the simulations using substate analyses. This was done by calculating the backbone dihedral angles of Met 20

Table 1
Differences in computed activation free energies, ⌬⌬G ‡ (kcal/mol), for WT ecDHFR and I14X mutants (X ‫؍‬ Val, Ala, Gly) relative to the WT with the Met 20

loop in the closed conformation at 25°C
The enzyme models were constructed with the Met 20 loop either closed or disordered. The corresponding free energy barriers are compared with relative experimental pre-steady state single-turnover data at 25°C. Calc., calculated; Exp., experimental.

Role of the Met 20 loop in the hydride transfer in E. coli DHFR
reported by Rod et al. (26). In the I14G* mutant, the Met 20 loop is highly disordered with a wide distribution of / angles. The impact of Met 20 on structure-The position of the GS stationary point of the PMF along the antisymmetric stretch coordinate for the WT enzyme in the closed conformation is approximately Ϫ1.8 Å, whereas for all other systems studied herein the value is greater, approximately Ϫ2.0 Å (Tables 2  and 3). There is, however, greater variability in the positions of the GS and TS for the mutant systems with a disordered Met 20 loop relative to the WT enzyme. For instance, the position of the TS ranges from Ϫ0.09 to Ϫ0.10 Å for the I14X systems (Table 2), whereas for the I14X* mutants the TS ranges from Ϫ0.05 to Ϫ0.16 Å (Table 3). Interestingly, the ensemble-averaged value of the DAD distance at the classical TS is nearly identical (2.63-2.66 Å) for all simulations. This invariability of the TS position with mutations is in agreement with earlier findings in our group (16,27) and many other groups (28 -34). However, the DADs in the GS show greater variability for the disordered structures than for the closed structures. The ensemble averaged D-H-A angle in the GS is rather similar for all systems with a closed Met 20 loop, ranging from 134 -145°, whereas this range is widened to 109 -149°in the disordered loop enzymes. However, ѯD-H-A at the TS is similar for all systems, ranging from 153 to 158°with similar standard deviations.
In accordance with MD simulations of the GS in Ref. 15, the interaction between the pterin ring and Asp 27 remains tight throughout the reaction for all systems. Similarly, the hydrogen bonds between the nicotinamide carboxamide moiety and the amide backbone of Ala 7 and Ile 14 remain intact for all systems except for I14G*. In this latter system, a different binding mode is observed in the crystal structure where the carboxamide group interacts with the backbone of Thr 123 instead of Ile 14 . Additionally, the WT enzyme with the closed Met 20 loop conformation has the tightest hydrogen bonds to the cofactor and substrate at the TS. The greatest difference in the interaction pattern between the WT/I14X and WT*/I14X* mutant systems is found in the contact between the Met 20 residue and the substrate-cofactor complex (Fig. 1). In the WT and I14X systems, the Met 20 loop packs against the reactive complex. This is reflected in the distances between the Met 20 sulfur atom and key nitrogen atoms in the NAD(H) and pterin moieties that range between 3.6 and 4.7 Å (Tables 2 and 3). In contrast, in the WT* and I14X* systems, the Met 20 loop is not packed against the pterin ring, and the range of Met 20 -nitrogen distances is 4.7-15.9 Å.
An additional interesting structural aspect that is greatly influenced by the conformation of the Met 20 loop is the active site hydration. It is well known that several conserved water molecules are located near the cofactor-substrate complex (6,(35)(36)(37)(38). Particular attention has been devoted to the water molecule likely responsible for the N5 protonation (38). Here we focus on a triad of waters that are located near the pterin ring O4 oxygen and Asp 27 (W1-3; Fig. 4) whose detailed structure has been revealed by a high-resolution WT ecDHFR structure (38). These water molecules stabilize the reduced pterin and Table 2 Ensemble-averaged structural properties of WT DHFR and I14X mutants with the Met 20

loop in a closed conformation at 25°C in the ground and transition states
Standard deviations are shown in parentheses. The ensemble-averaged data were obtained from 250-ps simulations of the GS and TS. The position of the TS along the reaction coordinate was determined from the position of the maximum PMF value. I14V  I14A  I14G   GS  TS  GS  TS  GS  TS  GS (40) found that this resonance-assisted charge-separated moiety (Scheme 2) is the most polarized region of the substrate-cofactor entity. In Fig. 5, we present the radial distribution functions (RDFs) around O4 of the pterin ring in the GS and at the TS for all enzyme forms studied herein. There is a striking difference in the RDFs for the closed and disordered Met 20 loop forms. In the closed enzyme, there are distinct peaks

Role of the Met 20 loop in the hydride transfer in E. coli DHFR
with the first solvation shell corresponding to a highly organized water molecule. The maximum of this peak is located at ϳ2.85 Å and corresponds to a water molecule that interacts with O4, Asp 27 , and Trp 22 (W1) and is highly conserved in DHFRs (Fig. 4). This water molecule is not exchanged throughout the simulations. A second peak in the RDFs is located at ϳ4.6 Å, and a third peak is found at ϳ7.5 Å, corresponding to the two water molecules W2 and W3 (see supplemental Fig.  S4). In the simulations with the disordered Met 20 loop, the RDFs are much broader, suggesting less hydration structure with considerably more waters penetrating the active site. However, even in the case of the open Met 20 loop simulations, a water molecule is directly interacting with O4. Inspection of the identity of the water molecule closest to O4 revealed that in the open Met 20 loop simulations this water molecule is exchangeable. Notably, in the GS of the I14G* simulations, the maximum peak is shifted to 2.95 Å, suggesting a more disordered state. For the WT enzyme with a closed Met 20 loop at the TS, the maximum peak of the first hydration shell is shifted to 2.75 Å, suggesting TS stabilization.
The impact of Met 20 on N5 pK a -We examined the pK a of substrate N5 atom in the ternary Michaelis complex. Our calculated pK a for the closed (WT) and disordered (WT*) state enzyme are 8.3 and 6.8, respectively, in excellent agreement with the values determined by Brooks and co-workers (14,41).The pK a obtained for the less ordered conformation (6.8) is similar to the experimentally reported pK a (ϳ6.5 (42) or 6.7 (43)), suggesting that protonation occurs in a conformational state with the Met 20 loop accessible to bulk water. In contrast, the pK a of N5 in the closed Met 20 conformation is observed to be less acidic (pK a ϭ 8.3), which is consistent with the recent experimental study by Liu et al. (43) These authors found DHF N5 to be partially protonated at pH values up to 11 based on kinetic studies (43). These kinetic studies showed that hydride transfer occurs even at such high pH values. Considering that in the current work we have established that hydride transfer occurs in the closed Met 20 loop state, we suggest that the elevated pK a value reported by Liu et al. (43) corresponds to the closed Met 20 loop state. The role of Met 20 in modulating the N5 pK a has been emphasized previously by Brooks and co-workers (14,41), and the current findings fully support their findings. We suggest that different Met 20 loop conformations bring about changes to the active site electrostatics that might affect N5 acidity. Specifically, these changes are expected to be correlated with the active site hydration as will be shown below.
Next, we analyzed the hydration near DHF in the closed and disordered Met 20 conformations as well as in bulk water. We plotted the radial distribution function between the substrate N5 and the water oxygens for the two folate forms, i.e. H 3 F ϩ and H 2 F (Fig. 6). The broad peaks near the H 3 F ϩ N5 (2.95 Å) in solution and WT* enzyme indicate specific hydration (Fig. 6A). In the closed WT complex, the Met 20 side chain packs against the pterin ring of the substrate (6,10,44,45), and no first solvation shell is observed. A small peak near 5.5 Å in both ternary complexes represents a water molecule stabilized by Asp 27 and H 3 F ϩ O4 interactions. In Fig. 6B, a sharp peak is observed at 2.85 Å, which shows that water tends to stabilize near N5 of H 2 F via hydrogen bond interactions in all the systems. The Met 20 loop flexibility has been suggested to be responsible for allowing water access to the substrate. Interestingly, in crystal structures of the ternary complex of ecDHFR, a water molecule is also observed near N5 of H 2 F (38). Seemingly, prior to N5 protonation, the packing of the Met 20 loop is somewhat loose, allowing access to water, whereas following N5 protonation the packing is tightened in preparation for the hydride transfer step, thereby excluding water molecules and hence raising the N5 pK a .

Role of the Met 20 loop in the hydride transfer in E. coli DHFR
The impact of Met 20 on dynamics-In Fig. 7, we present the root mean square fluctuations of the C␣ atoms as a function of amino acid residue number in the GS and TS for WT and I14X mutant enzyme forms in the closed and disordered Met 20 loop conformations. The main characteristic that is apparent in the root mean square fluctuation plots is the increased flexibility in the systems with the disordered Met 20 loop. In particular, the Met 20 loop (residues 9 -23) and F-G loop (residues 117-131) regions are considerably more flexible in the disordered than in the closed structures. This is likely due to changes in the contact between these two segments in the disordered structures. For instance, the contact between the side chain of Asp 122 in the F-G loop and the backbone amide of Glu 17 in the Met 20 loop is lost. We also note increased flexibility in the C-D loop (residues 64 -71) for all the disordered systems, including the WT* enzyme. These fluctuations are slightly greater in some of the mutants. Finally, we note that there is great flexibility in Gln 146 in the I14V* mutant in the disordered conformation. This greater fluctuation in the I14V* mutant is the result of the loss of an intraloop hydrogen bond between Asp 144 and Gln 146 during the simulation. The general observations discussed above are in agreement with earlier experimental and theoretical studies (11,46).

Discussion
A crucial question regarding the chemical steps in enzyme reactions is to what extent global protein motions might participate in the reaction coordinate. In the case of ecDHFR, it has been argued that the Met 20 loop participates in the reaction coordinate. In the current study, we asked an even simpler question: whether the Met 20 loop is preferably closed during hydride transfer or not and whether this influences the DADs. Seemingly, based on the free energy profiles obtained in this work, the Met 20 loop is indeed likely to be in the closed conformation in the WT enzyme with the Met 20 residue packed against the pterin ring. In the case of the WT enzyme, we simulated both the closed and disordered Met 20 loop configurations and obtained a slightly higher free energy barrier with the disordered Met 20 loop (⌬⌬G ‡ ϳ 1 kcal/mol). In the case of the I14X mutants (X ϭ Val, Ala, Gly), we obtained best agreement with pre-steady-state single-turnover rates using the more disordered structures (Table 1), although the results are within statistical error for X ϭ Val, Ala. A possible interpretation of these results is that the Met 20 loop in general and the Met 20 amino acid in particular are preferentially in the closed conformation during the hydride transfer step in WT ecDHFR but not necessarily so in the I14X mutants. That is, the I14X mutants can reach a reactive configuration similar to that of the WT enzyme but with a more disordered Met 20 loop. We note, however, that the effect of the Met 20 loop is small in agreement with experiments showing that the enzyme is active with significant parts of the Met 20 loop deleted (7). We did not observe any strong connection between the Met 20 loop conformation and the DADs.

Role of the Met 20 loop in the hydride transfer in E. coli DHFR
Our simulations verify earlier results of Brooks and co-workers (41) relating to the crucial role of the Met 20 loop in determining the substrate N5 pK a . The pK a value computed for the disordered Met 20 loop (ϳ6.8) is comparable with the experimentally determined value (ϳ6.5 (42) or 6.7 (43)), whereas the pK a computed for the closed conformation (8.3) is in agreement with the finding that hydride transfer occurs at basic pH values up to ϳ11 (43). The difference in pK a values for the two loop conformations is likely due to differences in active site hydration.
Additionally, the closed Met 20 loop conformation is also involved in ordering in the structure of several conserved water molecules in the vicinity of the pterin ring O4 position. It has been suggested previously that these water molecules are part of a proton relay (36,37) or ligand recognition (6,35,38); here we suggest that these water molecules also play a classical enzyme catalytic role, TS stabilization. Indeed, these tightly bound waters likely facilitate substrate polarization as suggested by Garcia-Viloca et al. (40) and stabilize charge separation (Scheme 2). This in turn enhances substrate binding via interaction with Asp 27 and serves to mitigate hydride charge transfer by screening the Asp 27 charge (Fig. 4). Experimentally, mutation of Trp 22 , which forms part of this hydration motif, to His or Phe results in a non-negligible reduction in the rate of the hydride transfer step as determined by stopped flow kinetics (47). In the case of other DHFRs, this could be a similar requirement, such as in human DHFR where, instead of the Met 20 residue, Leu 22 is packed against the substrate-cofactor complex and creates a similar water cavity.

Conclusions
A question concerning the catalytic cycle of ecDHFR that has garnered significant interest is whether the Met 20 loop motion is coupled to the hydride transfer step. In the current study, we examined whether the Met 20 loop is preferably closed during hydride transfer or not and whether this influences the DADs. To address this question, we studied the WT and three mutant forms (I14X; X ϭ Val, Ala, Gly) of the enzyme and compared the finding with published experimental observations. These enzyme forms were modeled with the Met 20 loop in its closed and disordered conformations. Using hybrid QM/MM free energy simulations, we obtained optimal agreement with experimental pre-steady-state kinetics modeling the WT in the closed Met 20 form and the mutants in the disordered Met 20 form. Indeed, we found that the disordered conformation of the Met 20 loop causes a consistent slight increase in the activation free energy barrier for the WT and mutant enzymes. We conclude that a closed Met 20 loop is likely preferential for catalysis as it aids in the preorganization of several conserved water molecules, likely influencing both ligand binding and TS stabilization. The DADs at the TS, however, remain rather similar in the closed and disordered conformations. We also found a significant correlation between the substrate N5 pK a and the Met 20 conformation. The open Met 20 loop conformation facilitates protonation prior to the hydride transfer. Closing of this loop reduces the acidity of the N5 proton, thereby stabilizing the H 3 F ϩ state in preparation for hydride transfer.

Materials and methods
The X-ray crystal structures of WT ecDHFR (Protein Data Bank codes 4RGC and 1RB2) (6,38) and mutant E. coli I14V, I14A, and I14G (Protein Data Bank codes 4QLG, 4QLE, and 4QLF, respectively) (16) with folate and the oxidized cofactor NADP ϩ were used to construct the initial configurations for the present study. The 4RGC structure corresponds to the closed ternary complex where the Met 20 residue is packed against folate (henceforth WT), whereas in 1RB2 the Met 20 loop is disordered (henceforth WT*; * ϭ disordered). The crystal structure of I14G was missing a part of the cofactor and part of a loop (residues 122 and 123). The missing cofactor coordinates were obtained by superimposing the I14V structure on the I14G structure followed by extraction of the missing coordinates for I14G from I14V. The missing residues were modeled with Discovery Studio software (Biovia, Inc.). These mutant enzyme forms that have the Met 20 loop in a disordered conformation are termed I14X* (X ϭ Val, Ala, Gly and * ϭ disordered). Additionally, a second set of mutant structures that are similar to the WT were constructed in silico by changing the Ile 14 residue in the original WT ecDHFR structure (Protein Data Bank code 4RGC) to Val, Ala, or Gly (termed I14X where X ϭ Val, Ala, Gly).
The setup and MD heating and equilibration simulations of all systems were carried out using procedures similar to those used for the WT and mutant ecDHFRs in previous studies (16,21,48). Briefly, the protonation states of all polar amino acid residue side chains were adjusted to pH 7, and the protonation states of the His residues (either neutral tautomeric forms or positively charged form) were determined based on the hydrogen bonding patterns of the local environment. The N5 position of the pterin ring of the dihydrofolate was protonated (43,49,50). The HBUILD facility in the program CHARMM was used to add hydrogen atoms (51,52). 29 sodium ions and 15 chlorine ions were added to all enzyme systems to neutralize the overall negative charge. This ionic concentration mimics experimental conditions (16) and effectively screens the charges in the system. The potential energy surface in the current study is described by a hybrid QM/MM Hamiltonian (53) where the QM region is treated by a modified AM1 semiempirical Hamiltonian (54) denoted AM1-specific reaction parameter (55). This Hamiltonian was designed to reproduce high-level calculations for an assortment of electronic and thermodynamic properties for reactions involving various nicotinamide and pterin derivatives (21). The QM region includes significant fragments of DHF and NADPH, which are proximal to the reaction center, whereas the MM region contains the remaining ligand atoms, the entire protein, water molecules, and salt. The water molecules were represented by the three-point charge TIP3P model (56). Hydrogen link atoms were placed across the bonds intersecting the QM/MM boundary. The QM/MM interactions were treated by electrostatic embedding. A detailed QM/MM partitioning scheme and a thorough description of the development of the AM1-specific reaction parameter Hamiltonian are provided elsewhere (21). In modeling the MM region, we used the all-atom CHARMM36 force field (57-60).

Role of the Met 20 loop in the hydride transfer in E. coli DHFR
Periodic boundary conditions were used to solvate the Michaelis complex using a pre-equilibrated cubic water box (ϳ65 ϫ ϳ65 ϫ ϳ65 Å) while treating long-range electrostatic interactions with the Ewald summation technique (64 ϫ 64 ϫ 64 fast Fourier transform grid, ϭ 0.340 Å Ϫ1 ) (61). All systems were fully minimized, heated up gradually to 298 K for 25 ps, and equilibrated at that temperature for 1 ns. All equilibrations were conducted in the isothermal-isobaric (NPT) ensemble at 1 atm, and the target temperature was controlled by the extended constant pressure/temperature method (62, 63) and a Hoover thermostat (64). The leapfrog integration scheme (65) was used to propagate the equations of motions, and the SHAKE algorithm (66) was applied to constrain all MM bonds involving hydrogen atoms, allowing a time step of 1 fs. During the initial stages of the equilibration, several nuclear Overhauser effect (NOE) restraints were imposed on key hydrogen bond interactions between the ligands and the surrounding residues as well as within the protein.
The umbrella sampling technique (67) was used to determine the classical-mechanical PMF for the hydride transfer reaction at 25°C (i.e. the free energy profile along a predefined reaction coordinate). The reaction coordinate was defined as the antisymmetric reactive stretch coordinate, asym . Specifically, asym is defined as the difference between the lengths of the breaking C4N-H and forming H-C6 bonds. The reaction coordinate was discretized and divided into 16 evenly spaced regions, or "windows," ranging from Ϫ2.00 to 1.5 Å. Each window was subjected to an appropriate harmonic restraint, which keeps asym in the desired region, and an umbrella potential (roughly the negative of the PMF) as a function of asym . The cumulative simulation time per window was 1 ns, and the statistics for the reaction coordinate were sorted into bins of width 0.01 Å. PMF profiles were computed using the weighted histogram analysis method (48,68). The statistical error was estimated using the bootstrapping algorithm (1000 steps) combined with the weighted histogram analysis method (69). The resulting error analysis of the free energy profiles obtained is provided in the supplemental Fig. S1. The effect of bin size on the free energy profiles was also investigated (supplemental Fig. S2). We also plotted the free energy profiles computed for different time intervals in the same figure.
To evaluate the pK a for H 2 F at the N5 position, we carried out free energy perturbation (FEP) simulations of substrate molecule in bulk water and in two ternary complexes of the WT ecDHFR (with closed and disordered Met 20 loops). We used single-topology FEP QM/MM simulations where a molecule is transformed from one chemical state to another via a series of intermediates (70,71). We used the -dynamics methodology where the potential energy function of system is defined as follows.
U͑͒ ϭ ͑1 Ϫ ͒U A ϩ U B (Eq. 1) The A and B states here correspond to the protonated substrate H 3 F ϩ and its unprotonated analog H 2 F, respectively, and is a coupling parameter. This approach entails annihilating the QM and classical force field terms for the proton attached to the DHF N5 position as a function of . In the present study, we chose ⌬ ϭ 0.1, and 10 windows were used to perturb ϭ 0 3 0.1 3 0.2 … 3 1.0. The cumulative simulation time per window was 1 ns. The double-wide sampling used here includes the FEP calculations in both forward and reverse directions. Therefore, for each sampling window (), the free energies for the forward and backward perturbations were obtained. We confirmed the convergence of these simulations by comparing the backward and forward perturbation free energies for each window. Finally, the protonation free energies for the full reaction were obtained by integrating all the sampling windows ().
The difference between the substrate protonation free energy in the enzyme and solution was used to compute the pK a shift for each complex relative to that in bulk water (72). To compute estimates for the experimental pK a values in the complexes, we added the corresponding computed ⌬pK a shift to the experimental pK a value in bulk water (2.6) (41,73).
Author contributions-D. T. M. designed and coordinated the study and wrote the paper. A. R. M. and A. V.-K. performed and analyzed the simulations. A. K. contributed to interpretation of the data and wrote the paper. All authors reviewed the results and approved the final version of the manuscript.