|Home | About | Journals | Submit | Contact Us | Français|
Conformational changes in enzymes are well recognized to play an important role in the organization of the reactive groups for efficient catalysis. This study reveals atomic and energetic details of the conformational change process that precedes the catalytic reaction of the enzyme dihydrofolate reductase. The computed free energy profile provides insights into the ligand binding mechanism and a quantitative estimate of barrier heights separating different conformational states along the pathway. Studies show that the ternary complex comprised of NADPH cofactor and substrate dihydrofolate undergoes transitions between a closed state and an occluded state via an intermediate “open” conformation. During these transitions the largest conformational change occurs in the Met20 loop of DHFR and is accompanied by the motion of the cofactor into and out of the binding pocket. When the cofactor is out of the binding pocket, the enzyme frequently samples open and occluded conformations with a small (≈ 5 kBT) free energy barrier between the two states. However, when the cofactor is in the binding pocket, the closed conformation is thermodynamically most favored. The determination of a profile characterizing the position-dependent diffusion of the Met20 loop allowed us to apply reaction rate theory and deduce the kinetics of loop motions based on the computed free energy landscape.
Dihydrofolate reductase (DHFR) is a biomedically important enzyme targeted for several anticancer and anti-bacterial therapies.1 Functionally, it catalyzes the reduction of dihydrofolate (DHF) to tetrahydrofolate (THF) with help of a cofactor nicotinamide adenine dinucleotide phosphate (NADPH). The product of this reaction, THF, is a precursor of cofactors required for purine, pyrimidine and amino acid synthesis.2 Because of its biological significance, DHFR has been subject of several theoretical3 and experimental4 investigations.
Based on the multiple crystal structures5 and kinetic data6 a catalytic pathway for DHFR has been deduced (Scheme 1). These studies suggest that during catalysis, DHFR cycles through five major intermediates: E:NADPH, E:NADPH:DHF, E:NADP+:THF, E:THF and E:NADPH:THF and undergoes significant conformational changes during this process. Analysis of the available crystal structures illustrates that the large-scale conformational change of DHFR is concentrated in its Met20 loop (Scheme 1). Thus depending upon the conformation of the Met20 loop, different DFTFR states along the catalytic pathway are characterized as open, closed or occluded. In the closed state the Met20 loop stacks against the nicotinamide ring while in the occluded state the Met20 loop sterically hinders the cofactor from binding in the active site. The conformation of the Met20 loop in turn seems to depend on the ligands bound in the substrate and cofactor binding sites.5 In the holoenzyme E:NADPH and the Michaelis Menten complex E:NADPH:DHF, the Met20 loop adopts the closed conformation. While in the other three product complexes along the DHFR catalytic pathway, E:NADP+, E:THF and E:NADPH:THF, the Met20 loop assumes an occluded conformation.
Experimental kinetic studies have provided the free energy differences between the thermodynamic end states of DHFR and the rates connecting them, but the free energy barrier heights are not known. Here, we have performed atomically detailed enhanced sampling simulations to determine the free energy profile and pathway corresponding to the conformational transition of the Met20 loop between the closed and occluded states of DHFR before the hydride transfer reaction (step 2, Scheme 1). Our results show that in the Michaelis Menten complex the closed state of the enzyme is thermodynamically favored over the occluded conformation by about 3.3 ± 1 kcal/mol. The highest free energy barrier corresponding to this conformational change is 5 ± 1 kcal/mol and the transition occurs via an intermediate ‘open’ conformation along the pathway. The low free energy barriers between different states along the pathway suggests that the ligand binding most likely occurs via selected-fit mechanism7 as recently proposed for other enzymes, e.g., adenylate kinase.8
The thermodynamic preference of the conformational states along the pathway obtained from simulations is in the excellent agreement with NMR studies.4 However, the computed free energy barrier is lower than that estimated using kinetic rates from NMR and conventional transition state theory with the Arrhenius type dependence of the rates. In order to resolve these differences, we determine the position-dependent diffusion coefficients of the Met20 loop and used Kramer’s reaction rate theory to calculate the rate of the Met20 loop transition in conjunction with the barrier height obtained from the free energy profile. The resulting rate of the Met20 loop transition is quite close to that obtained from NMR studies, confirming that different conformations of DHFR are indeed separated by the low free energy barriers. Overall, these results provide the missing dynamics picture between the available structural information and kinetic data corresponding to the conformational change process of the biologically important enzyme DHFR.
The free energy profile corresponding to the Met20 loop conformational change as a function of the ΔDrmsd order parameter is shown in Scheme 2. The closed state corresponds to the free energy minimum at ΔDrmsd = −3.5 Å. The occluded state corresponds to the ΔDrmsd = 2 Å. The free energy minimum corresponding to the occluded state is higher than that corresponding to the closed state by about 3.3 ± 1 kcal/mol. The nature of the ΔDrmsd progress variable is such that it allows structures to evolve away from the endpoint conformations simultaneously (supporting information, Scheme. S1) and helps identify possible stable intermediates along the pathway.9 As shown in Scheme 2, there is a stable intermediate conformation (ΔDrmsd ≈−l Å) of the Met20 loop along the pathway. This intermediate state has the characteristics of the crystallographically observed ‘open’ conformation of DHFR in which the Met20 loop faces away from the reactants resulting in the widening of the cofactor binding cleft.5 The free energy minimum corresponding to the open state is higher than that corresponding to the closed state by about 2.2 ± 1 kcal/mol. There are two free energy barriers along the Met20 loop conformational transition pathway. The barrier at ΔDrmsd = −1.4 Å, 4.5 ± 1 kcal/mol higher than the closed state, separates the closed state from the open state. Another, broad energy barrier in the ΔDrmsd range of 0.5–1.3 Å, 5 ± 1 kcal/mol higher than the closed state, separates the open state from the occluded state.
Transition between the closed and occluded states and vice-versa is accompanied by the movement of cofactor NADPH into and out the binding pocket. The thermodynamically favorable conformation of the Met20 loop corresponding to the position of the cofactor along the pathway can be traced from the two-dimensional free energy surface of the distance between the N5 atom of NADPH and C6 atom of DHF, involved in the hydride transfer reaction, with respect to the ΔDrmsd progress variable. As shown in Scheme 3, when the cofactor is out of the binding pocket (distance > 6 Å), two distinct stable conformations of the Met20 loop are observed. One of the conformers corresponds to the occluded state (ΔDrmsd ~ 2 Å) while the other conformer around ΔDrmsd ≈ −1.0 Å corresponds to an intermediate “open” conformation. The existence of these two conformations in the binary folate-bound form of DHFR was also suggested by high-pressure NMR studies.10 The cofactor most likely binds in this open state because of the wide active site cleft that provides access to the cofactor. Apparently, the open state is also the most frequently observed conformation crystallographically.5 However, when the cofactor enters the binding pocket, the closed conformation becomes thermodynamically most favored (Scheme 3). In the closed state, the distance between the atoms of the cofactor and DHF involved in the hydride transfer fluctuates in the range of 3–4 Å, which is conducive for the hydride transfer reaction that follows the conformational change step in DHFR’s catalytic cycle (Scheme 1).11 As shown in Scheme 3, there is only a small range of ΔDrmsd values (−3.5 Å– −2.8Å ) for which the cofactor and substrate are suitably placed for the hydride transfer. The Met20 loop conformations in this ΔDrmsd range can be classified as “tightly closed” loop states. Clearly, the Met20 loop can populate several energetically favorable conformational states (Scheme 2(a)). However, only certain conformations lead to overall active site geometry conducive for catalysis. Corroborating our findings, free energy perturbation simulations investigating the effect of the Met20 loop conformation on the substrate pKa concluded that it is only in the tightly closed loop state that the suitable pKa of 6.5 for the substrate protonation is attained.12
It has been suggested that protein dynamics is coupled to enzyme catalysis.3,13 In DHFR, mutation studies have suggested several side-chains of DHFR, in the Met20 loop, around the active site and even residues ≈20 Å away, play an important role in promoting catalysis.14 Notably, NMR studies of the E:folate:NADP+ complex revealed that the side-chains of residues Ilel4 and Ile94 populate both trans and +gauche rotamers about χ1 dihedral angle in solution.4 In the occluded and closed crystal structures only the +gauche conformation of these side chains is observed. Modeling suggests that in the trans rotameric state, the side-chain of these residues would overlap with the atoms of the cofactor and pterin rings, respectively. These observations led to the hypothesis that unfavorable interactions could be relieved by moving the cofactor and pterin rings of the substrate toward each other.4 In order to test this hypothesis we computed the free energy surfaces corresponding to the χ1 dihedral angle of the residues Ilel4 and Ile94 along the reaction coordinate, respectively (Scheme 4). As shown in Scheme 4, the side chains of both residues sample a minor trans (χ1 = 180°) rotamer state besides the dominant +gauche (χ1 = 60°) rotamer. The trans rotamer population of Ile94 lies in the ΔDrmsd range of 0.5 to −0.7 Å corresponding to the open state of the Met20 loop while that of Ilel4 lies in between −1 to −2.3 Å corresponding to the open and high energy conformations leading to the closed state (Scheme 2). Interestingly, the position of the trans population along the reaction coordinate coincides with a decrease in the hydride transfer distance (Scheme 3). Thus, the trans χ1 rotameric states of Ilel4 in fact steer the cofactor and folate toward each other, and likely facilitate catalysis. Additionally consistent with this suggestion, a mutagenesis study shows that the Ilel4 → Ala mutation in DHFR effects the hydride transfer rate by 20-fold.15 These observations reinforce the idea that the internal motion of the residue may be directly involved in rate-promoting vibrational modes and are coupled to the reaction coordinate. Interestingly, similar rate-promoting vibrations coupled to reaction coordinate were also suggested to facilitate hydride transfer reactions in the enzyme horse liver alcohol dehydrogenase.16
The rate of conformational transition from the closed to occluded state and vice-versa is available from NMR dispersion experiments.4 For the closed state, the rate of spontaneous fluctuations to the occluded state is 11 s−1 at 298 K. In the reverse direction, i.e., from the occluded to the closed state the rate is 440 s−1 at 298 K. According to conventional transition state theory (TST) the rate of reaction is related to the activation free energy as:17
where, k is rate constant, kB is the Boltzmann’s constant and T is the absolute temperature, h is Planck’s constant and ΔG† is the activation free energy. The expression kBT/h is a pre-factor whose numerical value near room temperature is ~ 6 × 1012 s−1. From simulations, the estimate for the rate of overall conformational change from the closed to occluded using standard transition state theory is available from Eq. 1. Using the frequency factor kBT /h and the highest free energy barrier, corresponding to the open to occluded Met20 loop transition, for which ΔG† = ~ 5 kcal/mol, the corresponding rate is 1.39 × 109 s−1. This estimated value is many orders of magnitude larger than the rate obtained from NMR.
Conversely, incorporating rate constants of the Met20 loop transition obtained from NMR4 into Eq. 1 yields an activation free energy barrier of kcal/mol from the closed to the occluded state and kcal/mol in the reverse direction. The free energy difference between the closed and occluded conformation is 2.2 kcal/mol , favoring the ternary closed state before the hydride transfer reaction. From simulations, the free energy difference between the closed and occluded states is 3.3 ± 1 kcal/mol with the free energy minimum at the closed state lower in energy with respect to the occluded state (Scheme 2). Thus, the free energy difference and the thermodynamic preference of closed and occluded states obtained from our simulations is in the excellent agreement with the NMR studies. However, the highest free energy barrier along the pathway is ≈5.0 kcal/mol, which is much lower than the transition state theory estimate of 16.0 kcal/mol from the closed to the occluded state.
In order to reconcile differences between the estimated rate and that obtained from NMR, we turn to Kramer’s reaction rate theory.18 Notably, the transition state theory estimate of the rate constant lacks dynamical information of the Met20 loop fluctuations. This missing information is incorporated into the rate estimates via the reaction rate theory. Kramer’s rate model assumes that the dynamics of the system can be represented as a diffusive process on a low-dimensional free energy surface. The representation of the complex dynamics of the system as diffusion along a well chosen reaction coordinate has been successfully applied to predict rates of protein folding19 and helix formation kinetics.20 Scheme 5, shows the position-dependent diffusion coefficients of the Met20 loop fluctuations obtained from umbrella sampling simulations. Interestingly, the highest value of the diffusion coefficient is found in the occluded state. Thus before hydride transfer, the conformations the Met20 loops adopts in the occluded state are restricted. There is faster diffusion from the occluded state toward the open state.
The intrinsic diffusion coefficient can used together with the free energy profile along the pathway to obtain an estimate of rate of loop transitions. The rate of reaction according to the Kramers-Smoluchowskii approximation is:18
where, and are the second derivative of the potential at the barrier top and in the bottom of the well, respectively, Dmax is the diffusion constant at the barrier top, ΔG† is the height of the free energy barrier, kB is Boltzmann’s constant and T is the absolute temperature (the effective mass is assumed to be unity for a fictitious particle diffusing on this free energy surface).
Assuming that the key ligand binding event occurs in the open conformation of the Met20 loop, we use the shape and height of the barrier separating the closed and open states to obtain a rate estimate (Scheme 2). Using Dmax values corresponding to ΔDrmsd ≈ −1.4 Å and a barrier height of 4.5 ± 1 kcal/mol into Eq. 2 yields a rate of transition from the closed state in the range of 23 × 104s−1 to 663 × 104s−1. This value is still larger than the NMR rates by about four orders of magnitude but smaller by the same magnitude from the crude transition state theory estimate. This is likely because our simulations employing an implicit solvent model are missing the granularity of the explicit solvent leading to smoother local effective energy surface. In general, explicit solvent has been shown to increase the roughness of the energy landscape by as much as 1.0 kcal/mol and slow down the dynamics of the system.21 Hamelberg et al. further suggested a rate correction factor of 104 at 300 K to account for the effect of explicit hydration on protein dynamics.21 On employing this correction factor, the rate of Met20 loop turns out to be 23 s−1 to 663 s−1. These values are in quite good agreement with NMR studies.4 More importantly, our results imply that the free energy barriers separating functionally important conformational states of the Met20 loop of DHFR are small and not large as suggested by including kinetic rates from NMR into the standard transition state theory. Thus, different conformational states along the pathway can be populated via thermal fluctuations.
Overall, the mechanistic picture emerging from the thermodynamic and kinetic analysis of Met20 loops fluctuations of DHFR is consistent with selected-fit mechanism of ligand binding.7 Specifically, when the cofactor is out of the binding pocket, the Met20 loop of the enzyme samples an ensemble of states comprising predominantly the occluded and open conformations and probably to some extent near closed conformations, and that are separated by small free energy barriers. As a cofactor approaches the binding pocket, the side-chains of the specific residues in the Met20 loop interact with the ring moieties of the cofactor and the substrate bringing them in a reaction competent state. Following, the Met20 loop stacks against the ring moiety of the cofactor resulting in the thermodynamically stable closed conformation.
An enhanced sampling method is used to obtain the free energy barrier heights of conformational change process preceding the catalysis in dihydrofolate reductase. The autocorrelation function of the order parameter for loop conformations is used to obtain position dependent diffusion coefficients. This information is used in Kramer’s reaction rate theory to obtain reaction rates. Our results show that different states along the Met20 loop conformational change pathway are not separated by large barriers as suggested from the naive application of transition state theory and rates from NMR. The observed low energy barriers between different conformations are consistent with the selected-fit model of ligand binding.
Two reference structures of DHFR were used to generate a path the between closed and occluded states. The Michaelis complex of DHFR in the closed conformation was modeled using the coordinates of the analogue ternary complex from PDB entry 1RX25 and the occluded conformation was modeled using coordinates from PDB entry 1RX6. In the occluded conformation, the missing coordinates of NADP were added from the complex of dihydrofolate reductase with 5,10-dideazatetrahydrofolate and nicotinamide adenine dinucleotide phosphate (PDB ID: 1RC4). In addition, the substrate analog 5,10-dideazatetrahydrofolate (ddH4F), was replaced with the coordinates of DHF in a configuration identical to that of the closed conformation. Hydrogen atoms were added to heavy atoms in both the closed and occluded conformations using the CHARMM HBUILD facility.22 The CHARMM program (c34al version)23 with the all-atom CHARMM22 force field was used for all calculations.24,25
In order to generate initial configurations for umbrella sampling, the Nudged-elastic band (NEB) method implemented in CHARMM26 was employed. Specifically, we used NEB to find a minimum energy path (MEP) between the closed and occluded states of DHFR. A total of 51 replicas including the end-point states were used to describe the conformational transition. Initial structures were generated by linear interpolation between the end-point states. Then the initial path was minimized with 500 steepest descent (SD) steps using the replica path method with a spring constant of 500,000 kcal/mol/Å2 without NEB force projection. The resulting path was minimized with NEB using adopted basis Newton-Raphson (ABNR) minimization for 2000 steps with the force constant of 10,000 kcal/mol/Å2. A final RMS off-path force of 0.002 kcal/mol/Å2 was reached. Further minimization did not change the gradient significantly. A distance-dependent dielectric coefficient was used to approximate solvent screening without including explicit water molecules during these initial stages of calculations. Electrostatic and van der Waals interactions were truncated at 16.0 Å.
Umbrella sampling MD simulations were performed to compute the free energy profiles along the conformational change path. The initial structures obtained from the NEB path were used as starting points for umbrella sampling MD simulations. The sampling was performed with a restraint Wj on the ΔDrmsd order parameter.9
where, ΔDmin is the value around which ΔDrmsd is restrained and Krmsd is the force constant. The ΔDrmsd order parameter is the difference in rmsd values of each structure from the reference reactant and product states. It is described as follows:
where, Xt is the instantaneous structure during simulation, Xclosed and Xoccluded are the two reference closed and occluded states of DHFR, respectively.
The 51 structures obtained from the NEB path optimization cover an rmsd range of 4.5 Å. These 51 structures formed 51 windows for umbrella sampling runs. In each window, the structure was subjected to restrained minimization using SD minimization for 1000 steps followed by ABNR minimization for 500 steps. During minimization a large harmonic restraint of 250 kcal/mol/Å2 was applied to protein heavy atoms in order to prevent the structures from drifting from the starting configurations. Then equilibration was performed with the restraint (equation 3) on protein heavy atoms, where the force constant was reduced from 210 kcal/mol/Å2 to 10 kcal/mol/Å2 over a period of 75 ps of constant temperature and volume MD simulation. Following the equilibration, production dynamics was performed and the structures were allowed to evolve with a weak one-dimensional ΔDrmsd restraint of 10 kcal/mol/Å2 on protein heavy atoms for 1.2 ns. Later, 7 more structures were added around ΔDrmsd between 0 to −1 Å in order to obtain uniform sampling. In this case structures were generated similarly using NEB. The restrained MD simulations were repeated thrice starting from the equilibrated structures with different set of initial velocities to ensure better statistics.
Equilibration and production dynamics was performed at constant temperature (300 K) using Langevin dynamics. The SHAKE algorithm was employed to constrain the bonds involving hydrogen atoms. This allowed the use of an intermediate time step of 0.0015 ps with the leapfrog integrator. A Langevin collision frequency of 10 ps−1 was chosen to couple the system to a 300 K heat bath. Electrostatic and van der Waals interactions were truncated at 16.0 Å with a switching function. The solvent effects were modeled through the GB/SA solvation scheme.27 Within this scheme, the electrostatic part of the solvation energy is computed from the generalized Born (GBMV) approximation. The non-polar contribution to solvation was added as ΔGnon–polar = γ·SASA, where SASA is solvent-acessible surface area, γ is set to 0.00542 kcal/mol·Å2. The missing CHARMM parameters of DHF and NADPH were borrowed from earlier work in our lab, where partial charges were determined by performing quantum mechanical calculations and missing atom types were assigned based on the analogy with other nucleic acids and/or amino-acids.12
The weighted histogram analysis method (WHAM)28 was used to obtain the potential of mean force along the one-dimensional ΔDrmsd reaction coordinate from the time series of the ΔDrmsd variable saved every 100 time-steps of 1.2 ns production dynamics. Two-dimensional free energy profiles along other additional degrees of freedom were generated using trajectories from the umbrella sampling MD simulations.9
Umbrella sampling simulations used to compute free energies of conformational transitions, as described above, were also used to compute diffusion coefficients along the ΔDrmsd reaction coordinate. The position-dependent diffusion coefficients (D(ΔDrmsd)) that describes the local dynamics of the Met20 loop on the free energy surface is given by following convenient expression,29
where, is the variance of the reaction coordinate ΔDrmsd in the harmonically biased simulations and is the characteristic time of its autocorrelation function
The autocorrelation functions were determined from the entire 1.2 ns restrained simulations with snapshots collected every 0.15 ps. The autocorrelation time at different ΔDrmsd values was determined by fitting a sum of two or three exponentials to its autocorrelation function.
This work was supported by the National Institutes of Health research grants RR012255 and GM048807.
Supporting Information Available
A media file depicting the conformational change of Met20 loop in DHFR and a figure showing the two-dimensional free energy surface along the reaction coordinate versus RMSD. This material is available free of charge via the Internet at http://pubs.acs.org.