Targeting a Rate-Promoting Vibration with an Allosteric Mediator in Lactate Dehydrogenase.

We present a new type of allosteric modulation in which a molecule bound outside the active site modifies the chemistry of an enzymatic reaction through rapid protein dynamics. As a test case for this type of allostery, we chose an enzyme with a well-characterized rate-promoting vibration, lactate dehydrogenase; identified a suitable small molecule for binding; and used transition path sampling to obtain ensembles of reactive trajectories. We found that the small molecule significantly affected the reaction by changing the position of the transition state and, through applying committor distribution analysis, showed that it removed the protein component from the reaction coordinate. The ability of a small-molecule to disrupt enzymatic reactions through alteration of subpicosecond protein motion opens the door for new experimental studies on protein motion coupled to enzymatic reactions and possibly the design of drugs to target these enzymes.

Allostery, the process by which biochemical processes are modulated by binding in a distant site, is complex and plays a major role in the maintenance of homeostasis in biological organisms. Historically, allostery was seen as a static effect, where an allosteric effector causes a structural change from one form, active or inactive, to the other.[1] One classic example is hemoglobin, where oxygen binding shifts the protein from the T to the R state, increasing oxygen affinity and leading to positive cooperativity.[2,3] Two phenomenological[4] models were developed to fit these observations of allostery, the Monod–Wyman–Changeux (MWC) model[5] and the Koshland–Nemethy–Filmer (KNF) model.[6] These static models for allostery were derived from stable structures for which X-ray crystal structures could easily be obtained, but the development of techniques such as NMR have provided dynamic structural information that allows for the extension of allostery to dynamic ensembles of states.[7]

The first theory of dynamic allostery was put forward by Cooper and Dryden[8] even before these dynamic structural methods were available and provides a statistical thermodynamic basis for entropically driven allostery. The authors state that cooperativity can originate from ligand binding which stiffens the vibrations of the protein, lowering the entropic barrier to the binding of subsequent ligands. Negative cooperativity would occur if the ligand loosens the vibrations of the protein. One example of a dynamically driven allosteric system is catabolite activator protein (CAP), which binds cyclic adenosine monophosphate (cAMP) with negative cooperativity.[9] CAP, which needs to bind two cAMP to allow for DNA binding and subsequent transcription, does not have an observable conformational change in the second cAMP binding site upon the first cAMP binding.[1,10] Because there is no change in the average structure, the negative cooperativity must arise from entropic effects.

Both static and dynamic allostery are changes in a binding affinity as a result of the binding of a ligand in another site. In this work, we further extend the concept of dynamic allostery, with small molecule binding that affects chemistry instead of substrate binding affinity. This is possible because some proteins exhibit a rate-promoting vibration (RPV), where subpicosecond protein motions play a role in the reaction coordinate.[11] Because bound molecules can alter the dynamics of the proteins to which they are bound, a small molecule in the right binding site could alter the RPV and disrupt the native reactive motions.

The concept of the RPV was first developed out of theoretical studies of proton exchange in crystals of benzoic acid, where symmetric coupling of the crystalline phonons to the reaction leads to rate enhancement.[12] This quantum theory was an extension of other theories of activated reaction rates in condensed phase, including Grote–Hynes[13] theory and others.[14,15] This theory was then applied to enzymes to identify protein vibrations that are coupled to the reaction by analyzing the spectral density of the protein.[16,17] These methods were crucial for identifying promoting vibrations in enzymes but do not provide atomistic scale information on the RPV. For that, transition path sampling is the optimal method.

Transition path sampling (TPS) is an unbiased Monte Carlo method which perturbatively samples the trajectory space of a rare transition between two stable, or metastable, regions of configuration space.[18−21] TPS provides rigorous thermodynamic ensembles, and calculation of the committor for trajectories in the ensemble can be used to obtain an ensemble of transition-state structures. The committor provides the probability of a certain configuration to reach the product region when dynamics is started with random velocities and the point where the probability is equal to one-half is defined as the transition state. Once an ensemble of transition states is obtained, committor distribution analysis can then be used to find the reaction coordinate perpendicular to the transition-state surface and the specific protein residues that participate, if any.[18,22] This technique has been used to identify RPVs in a number of proteins, including human dihydrofolate reductase,[23] purine nucleoside phosphorylase,[24] and lactate dehydrogenase.[25−29]

The RPV in human heart lactate dehydrogenase (LDH) in particular has been extensively studied. LDH catalyzes the conversion between pyruvate and lactate through reduction of a carbonyl and uses nicotinamide adenine dinucleotide (NADH) as a cofactor. The reaction occurs through a proton transfer from a histidine in the active site and a hydride transfer from NADH to the pyruvate. The reduction of pyruvate plays an important role in regenerating NAD+ in anaerobic conditions and serves as the terminal step of glycolysis in some bacteria and human cells.[2] To investigate RPV-mediated allostery, we chose to use LDH, because of its well-characterized RPV. We selected a binding site for the small molecule, docked a number of different molecules to this site, then performed TPS to generate two ensembles, one with a molecule bound and one without, to investigate allosteric effects. Computational details may be found in the Supporting Information.

To determine the ideal site for small-molecule binding, we identified an area which was close to the active-site, was adjacent to the RPV residues, and appeared to have a good accessibility from the solvent, so that an external molecule could easily bind. An image of the site location may be found in the Supporting Information. Using the crystal structure (PDB accession number 1I0Z),[30] we docked 158 molecules using Schrödinger Glide,[31] which were chosen based on their complementary structure to the chosen site. The molecule with the best docking score, 2-chloro-N-(3,5-dihydroxyphenyl)acetamide (CPA), shown in Figure , was chosen for further study. Systems with the bound molecule, called the CPA system, and without the bound molecule, called the control system, were prepared for TPS with explicit solvation, neutralization, minimization, heating, and equilibration.

Figure 1

Chemical structure of 2-chloro-N-(3,5-dihydroxyphenyl)acetamide, the small molecule bound to LDH as the possible allosteric effector.

When the structure of the two systems after equilibration was examined, CPA remained bound in the expected pocket and had somewhat altered the structure of the protein around it. This is shown in Figure . Especially notable is a shift in the position of the nicotinamide ring, such that the hydride donor–acceptor distance (DAD) in the CPA system was 0.28 Å greater than that in the control. A 0.5 ns dynamics run was performed to observe changes in the nonreactive dynamics of the protein. Root-mean-squared fluctuations (RMSF) for the short dynamics run were calculated and are shown in Figure . As is shown in the figure, most regions of the reactive monomer stiffen slightly upon binding of CPA. Major points of decreased flexibility include the regions around residues 41, 127, 152, and 178, which are all near the surface of the enzyme and distant from the active site. Residue 41 is at the opposite end of the same α-helix as a part of the RPV and the CPA binding site and may be directly affected by the bound molecule. The region of increased flexibility around residue 206 through 220 is a solvent accessible loop region. Arrows on the figure denote the three regions involved in the RPV in LDH, which consists of residues 31–33, 65–66, and 106.[22,26] The regions around the RPV are not significantly different in RMSF, which raises the question of whether the bound molecules is having an effect on the dynamics of the RPV at all. However, because a large RMSF is indicative of a flexible region with large motions and the RPV is a stiff, small motion, the RPV would likely not have a large RMSF. Additionally, the RMSF calculated is of equilibrium motions, which may be different from the protein motions during the reaction.

Figure 2

Selected residues from the equilibrated structure of the control (blue) and CPA (green) systems, aligned according to the quantum region. The bound molecule interacts with residues in the RPV, shifting thier position and slightly shifting the position of the nicotinamide ring.

Figure 3

Root mean squared fluctuations of the reactive monomer of LDH during 500 ps of dynamics. The control and CPA systems are shown. The black arrows denote regions previously identified as part of the RPV.

To observe the effect of CPA on the reaction, microcanonical TPS ensembles were generated with 200 trajectories in both the control and CPA ensembles. A selection of relevant average values for the control and CPA reactive trajectory ensembles are shown in Table . As can be seen from the table, the addition of CPA greatly increases the lag time between proton and hydride transfer while slightly decreasing both the proton and hydride DAD at the point of transfer. An increase in the lag time could negatively affect the rate of that step of the chemical reaction, but decreased DADs seem to indicate that CPA actually decreases the barrier to the independent particle transfers. Clearly these data are not enough to discern the net effect of CPA on the reaction.

Table 1

Table of Relevant Average Values for the Two TPS Ensembles, Including the Lag between the Proton and Hydride Transfer and the Average Proton and Hydride DAD at the Point of Proton and Hydride Transfer, Respectivelya

ensemble	lag time (fs)	proton DAD (Å)	hydride DAD (Å)
control	48.0 ± 30.4	2.56 ± 0.05	2.73 ± 0.03
CPA	124.7 ± 16.6	2.48 ± 0.05	2.66 ± 0.05

The point of transfer was defined as the point at which the particle is equidistant from the donor and acceptor.

Committor analysis was performed on every 20th trajectory in both TPS ensembles, with the main goal of collecting an ensemble of transition states for each ensemble of trajectories. Average values for the position of the transition state relative to proton and hydride transfer are shown in Table , and average values for the donor and acceptor values at the transition state are shown in Table . In the control system, the transition state of the reaction is concomitant with hydride transfer, similar to earlier studies on LDH.[22,29] In the CPA system, however, the transition state is decoupled from the hydride transfer; in fact, the transition state occurs closer to the proton transfer in this system. This could simply be a result of the increase in time lag between the two particle transfers or a change in the energetics of the system such that the hydride transfer step is no longer energetically dominant. The DAD in the transition-state ensembles are consistent with the transition-state locations, a small hydride DAD for the control system, but a significantly larger hydride DAD for the CPA system.

Table 2

Table of Average Position of Transition-States Relative to Particle Transfer from the Two TPS Ensemblesa

ensemble	PT-TSb (fs)	HT-TSc (fs)
control	45.7 ± 34.1	–3.6 ± 5.6
CPA	20.7 ± 33.1	–108.9 ± 28.8

Positive values are transfer before the transition state, and negative values are transfer after the transition state. As in Table , the point of particle transfer is when the particle is equidistant from the donor and acceptor.

Time between proton transfer and the transition state.

Time between hydride transfer and the transition state.

Table 3

Table of Average DAD at the Transition State from the TPS Ensembles with and without the Possible Allosteric Effector

ensemble	proton DAD (Å)	hydride DAD (Å)
control	2.68 ± 0.05	2.78 ± 0.07
CPA	2.64 ± 0.17	3.79 ± 0.24

To identify differences in protein contribution to the reaction coordinate, committor distribution analysis (CDA) was performed on four transition states from each ensemble. This method provides an approximate validity test for a guess of the reaction coordinate and was used to characterize the RPV in LDH.[22] Applying CDA to more than four transition states would be ideal, but CDA is computationally intensive, with approximately 300 000 quantum mechanical/molecular mechanical dynamics steps per transition state for each set of constraints. Distributions from CDA are shown in Figure . For the control system with the quantum region constrained, the distribution was peaked about 0.5 but had significant counts in the product and near the reactant well. This shows that constraining the quantum region is sufficient to constrain the majority of the reaction coordinate, but not all of it. When the protein is constrained along the RPV in addition to the quantum region constraints, the distribution is significantly improved, showing that protein motion along this coordinate plays a role in the reaction coordinate. In the CPA system, however, constraint of the quantum region yields a distribution that is peaked about 0.5, with only a small amount of distribution toward the product well. Additionally, constraining the protein along the direction of the reaction coordinate yielded a distribution that had a large amount of counts in the 0.7–0.8 region. This was surprising because, at least in theory, adding additional constraints to a set of constraints that fixes the reaction coordinate should not alter the distribution.

Figure 4

Histograms generated from committor distribution analysis, showing the control (a, c, and e) and CPA (b, d, and f) systems with three different sets of constraints. In panels a and b, the quantum region is constrained to the transition-state structure. In panels c and d, the protein is constrained along the RPV in addition to the quantum region with R106 constrained by the β-carbon. In panels e and f, the quantum region and RPV are constrained with R106 constrained by the ζ-carbon.

Examination of the individual constrained trajectories sheds light on the source of the difference in CPA committor distributions. When the RPV and quantum regions are constrained, the end of R106, which is adjacent to the pyruvate, moves away from the substrate. This is possible because the constraint on R106 is applied only to the β-carbon of the residue, which is distant from the guanidinium group that interacts with the pyruvate. When only the quantum region is constrained, the arginine also moves away, but to a significantly smaller degree. To test the influence of the arginine on the final committor distribution, we generated a committor distribution with the quantum region and the RPV residues other than R106 constrained in the same way as previously, but with R106 constrained to the substrate by the ζ-carbon, which is at the center of the guanidinium group. This distibution is also shown in Figure .

The distribution with these modified constraints in the CPA system recovered the distribution with just the quantum region. This shows that the distance from the guanidinium group to the pyruvate is in fact important, at least in the CPA system. We also applied the same modified RPV constraints to the control system, which resulted in a flat distribution. This indicates that the original RPV is a better reaction coordinate in the control system than the modified constraints. These distributions suggest that the small molecule changes the reaction coordinate to exclude the participation of subpicosecond protein motions that encompass large areas of the protein, but it does not exclude the local interactions of the pyruvate.

The increased hydride DAD in the equilibrated structure and removal of protein motion from the reaction coordinate suggest that the bound molecule increases the barrier to reaction. However, in seeming contradiction, the decreased DAD at proton and hydride transfer and the transition-state position farther from the hydride transfer suggest that the CPA decreases the barrier to reaction. A transition state significantly earlier than hydride transfer suggests that the barrier to hydride transfer no longer dominates the barrier to reaction. A crucial clue in understanding these data is that TPS can produce only ensembles of reactive trajectories. It does not provide information on the likihood of the reaction happening, only information on how it happens. Additionally, trajectories are populated according to their statistical weight within an ensemble, but between ensembles, the statistical weight of trajectories cannot be compared. This means that negative perturbations which do not completely eliminate the possibility of reaction will still be able to generate trajectory ensembles.

With this in mind, the decreased minimum DADs at transfer may not be due to an increase in reaction probability, but are likely due to a decrease in reaction probability eliminating trajectories with a longer distance. Here the donor and acceptor move on their own without mediation of the protein. Additionally, the transition-state shift could be due to a lower barrier to hydride transfer or other effects, such as an increased barrier to proton transfer or the increased time between particle transfers. While a lower barrier to hydride transfer lowers the total overall barrier to reaction, an increase in barrier to proton transfer and an increased time lag increase the barrier to reaction. Looking at all the data as a whole, it appears that the effect of the CPA on the reaction is entirely negative; the equilibrium structure is disrupted so that it is in a less reactive conformation, the reaction is forced to react at a smaller distance because of the removal of protein contribution, and the transition state is shifted with an increased time lag. This is the first representative of RPV-mediated allostery, where a small molecule disrupts the RPV of an enzyme to disrupt the chemical reaction.

Further study is required to gain more insight into allostery mediated through subpicosecond protein motions. Experimental studies would be instructive, though the molecule we chose, because this was intended as a proof of concept for RPV-mediated allostery, may have experimental issues such as solubility that were not considered. One informative experiment would be to measure the kinetic isotope effect on the hydride in H2O and in D2O. Determining whether the raised barrier for proton transfer in D2O suppressed the observed kinetic isotope effect for the hydride transfer would shed light on the change in position of the transition state in the CPA system. Other computational studies could be directed toward identification of the barrier to reaction with and without a small molecule. Unfortunately, typical barrier methods based on transition-state theory average over the protein motion and would thus fail to capture all of the protein dynamics in the reaction.[32] Barrier methods within the TPS framework have been used to observe the barrier for the transfer of specific particles but have not been extended to calculate a barrier for a full reaction.[33] Rate calculation methods have also been developed using TPS, but such methods are computationally expensive.[34]