Role of substrate dynamics in protein prenylation reactions.

CONSPECTUS: The role dynamics plays in proteins is of intense contemporary interest. Fundamental insights into how dynamics affects reactivity and product distributions will facilitate the design of novel catalysts that can produce high quality compounds that can be employed, for example, as fuels and life saving drugs. We have used molecular dynamics (MD) methods and combined quantum mechanical/molecular mechanical (QM/MM) methods to study a series of proteins either whose substrates are too far away from the catalytic center or whose experimentally resolved substrate binding modes cannot explain the observed product distribution. In particular, we describe studies of farnesyl transferase (FTase) where the farnesyl pyrophosphate (FPP) substrate is ∼8 Å from the zinc-bound peptide in the active site of FTase. Using MD and QM/MM studies, we explain how the FPP substrate spans the gulf between it and the active site, and we have elucidated the nature of the transition state (TS) and offered an alternate explanation of experimentally observed kinetic isotope effects (KIEs). Our second story focuses on the nature of substrate dynamics in the aromatic prenyltransferase (APTase) protein NphB and how substrate dynamics affects the observed product distribution. Through the examples chosen we show the power of MD and QM/MM methods to provide unique insights into how protein substrate dynamics affects catalytic efficiency. We also illustrate how complex these reactions are and highlight the challenges faced when attempting to design de novo catalysts. While the methods used in our previous studies provided useful insights, several clear challenges still remain. In particular, we have utilized a semiempirical QM model (self-consistent charge density functional tight binding, SCC-DFTB) in our QM/MM studies since the problems we were addressing required extensive sampling. For the problems illustrated, this approach performed admirably (we estimate for these systems an uncertainty of ∼2 kcal/mol), but it is still a semiempirical model, and studies of this type would benefit greatly from more accurate ab initio or DFT models. However, the challenge with these methods is to reach the level of sampling needed to study systems where large conformational changes happen in the many nanoseconds to microsecond time regimes. Hence, how to couple expensive and accurate QM methods with sophisticated sampling algorithms is an important future challenge especially when large-scale studies of catalyst design become of interest. The use of MD and QM/MM models to elucidate enzyme catalytic pathways and to design novel catalytic agents is in its infancy but shows tremendous promise. While this Account summarizes where we have been, we also discuss briefly future directions that improve our fundamental ability to understand enzyme catalysis.

Introduction

The relationship between a protein’s structure and its function has been the subject of intense interest for many years.[1] From the earliest crystal structures, mechanistic hypotheses could be formulated and then tested using a range of biochemical techniques. In the past decade the linkage between protein dynamics and catalysis has been of intense interest and controversy.[2−4] The notion here is that fluctuations in the protein “matrix” can notably affect catalytic rate. With the advent of site-directed mutagenesis, biochemists and biophysicists were able to further test the role of individual amino acids on the function of a protein.[5] While in aggregate all of these efforts enhanced our understanding of enzyme catalysis, our ability to rationally design a protein to be a highly efficient catalyst remains a significant contemporary challenge.[6−8]

Herein we describe two systems involving prenylation of various substrates. The first system is farnesyltransferase (FTase), which catalyzes the attachment of farnesylpyrophosphatase (FPP) to a cysteine in a conserved sequence Ca1a2X at or near the C-terminus of a protein.[9] While the reaction is facilitated by an active site Zn(II) ion,[10] the reaction rate is enhanced 700-fold via the presence of Mg(II),[11] whose binding site was initially unknown. The ground breaking structural work of the Beese lab provided a detailed mechanistic insight into several structural waypoints along the reaction pathway.[10] A remarkable feature of the experimental structures was the presence of FPP bound in the active site ∼8 Å from the cysteine bound to the active site zinc ion.[10] How this gap is overcome and how much it contributes to the overall barrier was unknown at the time. Through molecular dynamics (MD) simulations and combined quantum mechanics/molecular mechanics (QM/MM) simulations, we identified the putative binding site for Mg(II)[12] and the mechanism by which the active gap is overcome.[13] Finally, through detailed QM/MM simulations, we further rationalized the detailed kinetic isotope effect (KIE) from the Fierke lab by providing insights into the structure of the transition state (TS) for the rate determining prenylation step.[14]

The second system that we describe is NphB, which is an aromatic prenyltransferase (APTase) that has a mixed product distribution. Experimentally, the geranyldiphosphate (GPP) prenylation of 1,6-dihyroxynaphthalene (1,6-DHN) by NphB lead to three products with the first two being observed in a 10:1 (5-geranyl-1,6-DHN and 2-geranyl-1,6-DHN) product distribution and the third having a much smaller kcat than observed for the other two (4-geranyl-1,6-DHN).[15−17] Concomitant X-ray crystal structures of NphB with 1,6-DHN bound in the active site could only explain the major product.[15−17] The question of how the minor products arise and how would it be possible to favor these product over the major product through modifications of the active site initiated our interest in studying this system. Through MD[18] and QM/MM[19] simulations, we were able to demonstrate substrate dynamics involving multiple free energy wells that lead to the altered product distribution.

Recently we have also reported on the reactions catalyzed by CloQ[20] and FTmPT1[21] ATPases, that have interesting product distributions. For the latter it was observed that one simple point mutation (glycine to threonine) yielded a radically different product than that observed for the native protein.[22] These are interesting problems in and of themselves, but for the sake of this Account, we will focus on FTase and NphB since these are more complete stories.

Methods

To study substrate dynamics, we perform MD simulations[23] using an isobaric and isothermal ensemble (NPT) in explicit water while employing the particle-mesh Ewald method[24] for long-range electrostatics. In classical MD simulations, zinc ions were modeled using the popular “bonded model” approach that has been developed by us.[25−31] To map out the nature of the active site dynamics we also utilize potential of mean force (PMF) simulations,[23] where requisite degrees of freedom are held fixed (for example, a bond distance) while the remaining degrees of freedom are allowed to move. The PMF approach allows one to obtain the free energy cost for the restrained degree of freedom in both classical MD simulations and QM/MM studies of reaction barriers.

The combined QM/MM[32] method allows for the use of QM in the region of interest, which in our case is the active site and substrates that undergo reaction, while the remaining degrees of freedom utilize the less expensive classical potential like the AMBER force field[33,34] utilized here. The QM/MM approach[32] has found widespread applicability and its development was recognized by the 2013 Nobel Prize in Chemistry. The QM method used in our studies was the self-consistent charge density functional tight binding (SCC-DFTB) method,[35,36] which is a semiempirical model that combines good accuracy with relatively fast computational performance.[36] In our studies, we find that the SCC-DFTB method reproduces prenylation processes with excellent accuracy as judged by our ability to match experimental free energy barrier heights.[14,19] Technical details of the methods employed in our studies can be found in the extant literature and will not be described fully in this Account.[14,19]

Applications of MD and QM/MM to Prenylation Reactions

We have used MD and QM/MM methods to explore the structure, function, and dynamics of prenyltransferase reactions catalyzed by FTase,[12−14] NphB,[18,19] CloQ,[20] and FtmPT1.[21] Below we briefly highlight our efforts to elucidate the role of substrate dynamics and its effect on the FTase and NphB catalyzed reactions.

FTase

FTase and geranyltransferase (GGTase) are zinc metalloenzymes that catalyze the attachment of a farnesyl (a 15-carbon isoprenoid) or geranylgeranyl (a 20-carbon isoprenoid) group provided by farnesyldiphosphate (FPP) or geranylgeranyldiphosphate (GGPP) to a cysteine residue at or near the C-terminus of protein acceptors (see Figure 1).[9] This posttranslational modification is important for many GTP-binding switch proteins on receptor tyrosine kinase (RTK) signal transduction pathways and many proteins that are downstream, facilitating their attachment and localization to the inner side of the plasma membrane.[37−42] These proteins function as molecular switches, regulating cell proliferation and differentiation and cell survival and modulating cellular metabolism. Their malfunction is often associated with uncontrollable cell growth, which may lead to tumor and cancer formation. The catalytic mechanism for FTase is provided in Scheme 1. The binding of FPP occurs first (1 → 2) and is followed by the binding of the protein or peptide target[43,44] (2 → 3; see Scheme 1). The crystal structure of the FTase reactant ternary complex (protein data bank (PDB) code 1QBQ) finds a striking 7.4 Å gap between the C1 in the FPP substrate and the Sγ of the zinc-bound cysteine of the peptide target (see Figure 2 and Scheme 1), occupied by several solvent molecules. The farnesyl group of the native FPP substrate binds similarly in the FPP-bound FTase binary complex (PDB code 1FT2), as does HFP in the ternary complex (PDB code 1QBQ); hence, this gap is thought to be present in the reactive FTase ternary complex as well. Comparing the conformation of the farnesyl group in several FTase complexes with FPP or FPP analogs bound to the product complex (PDB code 1KZP) gave rise to a popular explanation of how the 7.3 Å gap is overcome. This proposal hypothesizes that a rotation of a bond between the first and second isoprene groups diminishes the gap between the C1 of FPP and the attacking Sγ and is supported by mutagenesis studies.[11]

Figure 1

Protein prenylation reactions catalyzed by FTase and GGTase.

Scheme 1

A Catalytic Mechanistic Hypothesis for FTase

Figure 2

Striking gap between the two reacting centers is demonstrated with a transparent “bond”. Also shown are the crystal water molecules found within 5 Å of HFP, the FPP analog. The zinc ion is shown as a yellow sphere, while oxygen atoms from water molecules are shown as red spheres.

The displacement of diphosphate from FPP by the zinc-bound thiolate (4 → 5) has been predicted to be SN1-like[45] as well as SN2-like.[46,47] Fierke and co-workers[43] predicted an associative/dissociative mechanism, which was also computationally observed by Klein et al.[48] Besides the role of zinc in the farnesylation reaction, Mg2+ ions play a supporting role in the reaction. The binding of Mg2+ ions at the millimolar level increases the reaction rate by 700-fold,[43] even though Mg2+ is not required for FTase function. The role of Mg2+ in enhancing the reaction rate is not clearly understood because there is a lack of structural information on how Mg2+ binds. Mutagenesis studies[11] suggest that magnesium is bound between the two phosphate groups of FPP and possibly interacts with D352β, stabilizing the leaving diphosphate group during farnesyl transfer. The overall rate of farnesylation is controlled by product release (5 → 6; see Scheme 1), which may be assisted by the binding of another substrate molecule. In mammals, the rate constant for product release is 0.05 s–1, while the rate constant of farnesylation is 17.0 s–1.

In the following, we summarize our studies aimed at (a) elucidating the binding position for Mg2 in the active site pocket of FTase, (b) understanding how the gap between FPP and the zinc-bound cysteine is overcome via conformational changes, and (c) the nature of the FTase prenylation transition state.

Mg2+ Binding to FTase

We have performed studies to help understand the role that magnesium plays in the catalysis of FTase.[12] We first modeled the process of Mg2+-binding in the active site of FTase using information from experimental studies of FTase and similar systems in order to generate starting configurations for our studies.[11] MD simulations were carried out to explore the validity of these assumptions. Through multi-nanosecond classical MD simulations, two types of Mg2+ binding positions were determined (see Figure 3): (1) the divalent metal ion interacts with two negatively charged oxygen atoms of the α-phosphate of the FPP3–, the carboxylate group of Asp147α, and three crystal water molecules; (2) the magnesium ion interacts with two negatively charged oxygen atoms from the α- and β-phosphate of FPP3–, the carboxylate group of Asp352β, and three exchanging water molecules. The second motif is in agreement with experiments performed by the Fierke group that find Asp352β to be important for the catalyic process.[11] Furthermore, in this motif, the distance between two reacting centers, the C1 atoms on FPP3– and the Cβ atom on acetyl-CVIM (5.2 Å), is shorter compared with the case when Mg2+ is not present (7.2 Å). Overall, our calculations suggest that Mg2+ ions may enhance the reaction rate in part by reducing the gap between the reacting partners.

Figure 3

Mg2+ binding motifs identified from MD simualtions. (a) Type 1, in which Mg2+ interacts with Asp147α, three water molecules, and both oxygen atoms on α-phosphate of FPP3–; (b) type 2, in which Mg2+ interacts with three water molecules, Asp352β, and two oxygen atoms from α- and β-phosphate of FPP3–. Zn2+ and Mg2+ are shown as green and silver spheres, respectively.

Overcoming the Gap

Prior to studying the chemical step in FTase, we aimed to understand the nature of the conformational change that brought FPP in close proximity to the zinc-bound thiolate. Indeed, isotope effect experiments suggested that this may be the rate-limiting step for some peptides.[43,49] The charged state of FPP proved to be a further complicating factor, with FPP existing as FPP3– or in its protonated form FPP2–. Extensive MD simulations in the absence of Mg2+ found that while the protonated form, FPP2–, may undergo the critical conformational change, the more negatively charged form, FPP3–, remains locked in place owing to the presence of positively charged residues in the active site pocket of FTase.[50] This effect can be seen in Figure 4, where PMF simulations bringing FPP closer to the zinc-bound thiolate are presented for both FPP3– and its protonated form for the CVIM peptide. From the solid line in Figure 4, we see the free energy cost for the conformational transition rise from a few kilocalories per mole in the FPP2–/CVIM case to becoming prohibitive in the FPP3–/CVIM case. In the presence of Mg2+ ions, FPP3– can also undergo this conformational change. Other theoretical groups have also explored the issues relating to bridging this gap.[51−53]

Figure 4

Free energy profiles of the conformational activation step in FTase ternary complexes computed as a function of the distance between the center of mass of atoms C1, C2, and O1 of FPP3– or FPP2– and Sγ of Cys2 of the target peptides. The full PMF curve of FTase/FPP3–/CVIM is truncated at 3.0 kcal/mol to fit into the plot of the other PMFs.

Having developed an understanding of the starting state in the absence of Mg2+, we next performed PMF simulations to evaluate the free energy for the conformational change associated with bringing the reactive centers into contact with one another for a number of scenarios. Our calculations found that for FPP2–/CVIM and the FPP2–/CVLM cases, the conformational transition cost only a few kilocalories per mole. In these PMF calculations, we observed that Y300β acted as the donor of the hydrogen bond with the β-diphosphate of FPP throughout the entire reaction profile. We hypothesized that the catalytic power of Y300β (∼500-fold decrease in kchem upon mutation) is not associated with the FPP activation step. To investigate this, we computed the corresponding free energy profile for the conformational activation of the Y300Fβ mutant FTase and CVIM (the dotted line in Figure 4). In support of our hypothesis, we found that eliminating the possibility of the β-diphosphate–Y300β hydrogen bond does not significantly alter the free energy of the conformational transition. We find that the mutation shifts both the intermediate state and the resting state farther away from the zinc-bound cysteine of the peptide. These results indicate that the measured reduction in the reaction rate upon this mutation is likely from the chemical reaction step in which Y300β may participate in the stabilization of the leaving diphosphate group. Overall, our calculations found the conformational change to have an activation barrier lower than that of the overall experimental rate constant, making it likely that the chemical step would be the rate-limiting step in most instances.

Structure of the Transition State to Prenylation

Following our analysis of the conformational change, we next examined the catalytic mechanism of the chemical step for FTase catalysis in great detail by performing QM/MM PMF calculations.[14] In calculations performed on the FTase/FPP2– system, the QM region consisted of the CVIM peptide chain, FPP2–, Zn2+, and the side chains of the four Zn2+ ligands. All other residues and water molecules were included in the MM region. We carried out a QM/MM PMF simulation from a C–S distance of 8.0–1.8 Å that mapped out the chemical step as well as the conformational transition.[13] The free energy profile accompanied by snapshots of the transition state and intermediates from our calculation is shown in Figure 5. The calculated barrier height of 20.6 kcal/mol is in excellent agreement with the experimental values for the GCVLS (20.0 kcal/mol) and TKCVIF (21.1 kcal/mol) peptides.[54,55] The QM/MM calculated conformational transition part of the profile is in excellent agreement with our previous MM based work[50] and revealed a shallow intermediate at 3.8 Å, representing a prereactive conformation. Furthermore, similar QM/MM PMF simulations performed on the mutant Yβ300F found that the mutation destabilizes the transition state by 1.8 kcal/mol, in good agreement with the experimentally determined value of 2.6 kcal/mol.[54,55]

Figure 5

FTase reaction profile including the conformational step (8.0–4.0 Å) and the chemical step (4.0–1.8 Å). The X-axis is the C–S distance (Å) and the Y-axis is free energy (kcal/mol). Snapshots of the intermediates and the TS are shown as insets.

In a collaborative effort with the Fierke group,[14] we examined the differing 3H α-secondary kinetic isotope effects for the CVIM and CVLS peptides. These two substrate peptides of protein FTase have been suggested to have two different rate-determining transition states (TSs) in the chemical step of the enzyme.[49] Based on 3H α-secondary kinetic isotope effect measurements, the former was proposed to have a rate-limiting SN2-like TS with dissociative characteristics, while due to the absence of an isotope effect, the latter was proposed to have a rate-limiting peptide conformational change.[49] Using PMF QM/MM simulations, we observed the experimentally proposed TS for CVIM but found that CVLS has a symmetric SN2 TS, which yielded a zero isotope effect, which is consistent with the experimentally observed 3H α-secondary kinetic isotope effects.[14] As such, our simualtions helped explain the mechanistic dichotomy as arising from changes in the TS structure.

Our FTase studies afford an excellent example of where computational efforts can supplement and enhance experimental insights. Over the course of these studies, we have expanded our understanding of the FTase catalyzed reaction by elucidating the role of the conformational change in the reaction mechanism. Our calculations have further aided in the interpretation of the experimentally observed kinetic isotope effects by demonstrating the lack of a rate-limiting step associated with the conformational change and further via the observation of a symmetric TS that yielded a zero isotope effect as determined experimentally.

NphB

NphB is an APTase that catalyzes the attachment of a 10-carbon geranyl group to aromatic substrates.[15−17] Importantly, NphB exhibits rich substrate selectivity and product regioselectivity as summarized for 1,6-dihydroxnaphthalane (1,6-DHN) in Figure 6.[15−17] We have performed a systematic computational study in order to understand several key questions related to the NphB catalyzed geranylation reaction.[18,19] In these studies, we first aimed to understand how three different products were formed over the course of the reaction, though the crystal structure of 1,6-DHN bound suggested a pathway for only one product (5-geranyl-1,6-DHN). Second, we attempted to test the hypothesis that the reaction mechanism of the prenylation step is a SN1 type dissociative mechanism with a weakly stable carbocation intermediate.[15−17] Finally, given the numbers of aromatic groups in the active site region, we wanted to determine the role of these residues in stabilizing potential carbocation intermediates along the reaction coordinate.

Figure 6

Schematic representation of geranylation catalyzed by NphB complexed with GPP and 1,6-DHN. From left to right, the product distribution is 10:1:trace.

To better understand the dynamics of 1,6-DHN in the active site of NphB, we performed extensive MD and MD PMF studies that calculated the free energy surface of this ligand ensconced in the binding pocket. The X-ray crystal structure (PDB code 1ZB6) placed the carbon atom that lead to the major product in closest proximity to the C1 carbon of GPP. The calculated free energy surface of 1,6-DHN bound in the pocket is shown in Figure 7. In this figure, we have not used the canonical carbon atom labels, but it can be readily oriented to Figure 6 by mapping C2 to canonical C5 and C9 to canonical C2. Two minima were identified from our calculations: the first minimum at the lower left of the right-hand panel in Figure 7 readily allows for the formation of the major product 5-geranyl-1,6-DHN. The second minimum is ∼2.3 kcal/mol higher in energy than the global minimum and is the starting point for the formation of the two minor products (e.g., 2-geranyl-1,6-DHN and 4-geranyl-1,6-DHN see Figure 6).

Figure 7

(a) The crystal conformation of the substrates (Mg2+, GPP, and 1,6-DHN) in 1ZB6 and (b) the computed 2D free energy contour plot (kcal/mol) using distances between C1 in GPP and C2 (canonical C5) and C9 (canonical C2) in 1,6-DHN. C1, C2 (canonical C5), and C9 (canonical C2) are represented as orange spheres.

The product distribution for a wide-range of products observed for the geranylation of 1,6-DHN catalyzed by NphB complexed with geranyldiphosphate (Figure 6) is partly due to the free energy preferences of the substrate binding states that favor one reaction channel over another.[15,16] Figure 8 shows the free energy profiles calculated from PMF simulation results using the SCC-DFTB QM/MM approach available in AMBER.[35] The excellent agreement between the computed and experimentally observed activation free energy values demonstrates that SCC-DFTB performs satisfactorily for this class of reactions.[19] The observed difference in the rates of product formation from 5- and 2-prenylation arises from the differing orientations of the aromatic substrate in the 1,6-DHN resting state (compare red vs blue starting point in Figure 8). While 4-prenylation shares the same resting state with 5-prenylation, the lower free energy barrier for carbocation formation makes the latter reaction more facile (blue vs. green curve in Figure 8, also see Figure 7), providing a rationale for why 4-geranyl-1,6-DHN is only found in trace amounts relative to 2 and 5 prenylation. Fianlly, the high free energy barrier associated with 7-prenylation is caused by the unfavorable orientation of 1,6-DHN in the active site pocket.

Figure 8

Free energy profiles for prenylation at four sites in 1,6-DHN.

Thus, the energy difference in the substrate binding position that favors one reaction channel over the other can explain the observed minor and major product distribution of 10:1. The MD simulations described above gave insights into the dynamics of the substrate within the binding pocket and afforded an explanation of the observed product distribution. QM/MM reaction path studies expanded upon the classical modeling of enzyme conformational dynamics and helped explore the reaction channel preferences.

A novel π-chamber composed of Tyr121, Tyr216, and the substrate 1,6-DHN was found to be important in stabilizing the carbocation intermediate (see Figure 9). Additionally, the π-chamber served to protect the intermediate by sequestering it away from water molecules. Our QM calculations (M06-2X/6031G** with basis-set superposition error (BSSE) corrections)[56−58] find that for the preferred 2- and 5-prenylation channels, the π-chamber stabilizes the geranyl carbocation by −20.6 and −13.4 kcal/mol, respectively. For the trace product arising from 4-prenylation and for the unobserved 7-prenylation pathway, the stabilization is significantly less (−6.8 and +1.0 kcal/mol). As such, our calculations suggest that the π-chamber serves a dual function: it protects the carbocation from water and selectively stabilizes the forming carbocation for the favored reaction channels.

Figure 9

π-Chamber in the NphB binding pocket consisting of Tyr121, Tyr216, and 1,6-DHN found at the resting state (left) and intermediate state (right) of 5-prenylation.

A facile water mediated proton transfer facilitates the loss of hydrogen at the prenylation site to form the final prenylated product in all cases (barriers are all below 10 kcal/mol) making carbocation formation the rate-limiting step. Interestingly, the same crystallographically observed water molecule was found to be responsible for proton loss in all three experimentally identified products (5-, 2-, and 4-prenylation). We find that after proton transfer, the relaxation of the final product from a sp3 carbon center to a sp2 center triggers a “spring-loaded” product release mechanism that pushes the final product out of the binding pocket toward the edge of the active site. The hydrogen bond interactions between the two-hydroxyl groups of the aromatic product and the side chains of Ser214 and Tyr288 help “steer” the movement of the product. In addition, mutagenesis studies[15−17] identify these residues as being responsible for the observed regioselectivity, particularly for 2-prenylation.

Our observations provided valuable insights into NphB chemistry, offering an opportunity to better engineer the active site and to control the reactivity in order to obtain high yields of the desired products. We find that substrate stabilization plays the key role in the case studied, but stabilization of the intermediate in the π-chamber could also play a role in mutant systems. Nonetheless, NphB is an interesting platform for future catalyst design work given the simplicity of the reaction and its level of experimental and theoretical characterization. Furthermore, the SN1 reaction mechanism observed for NphB differs from the prenylation reaction found in, for example, the farnesyltransferases (see FTase above), which proceeds via an SN2-like reaction pathway.[14] Many unanswered questions, such as the product release dynamics via the spring-loaded release mechanism, will have to be addressed to make the design of novel catalysts on the NphB platform a reality.

The Future of MD and QM/MM Methods in Catalyst Design

MD and QM/MM methods have had a major impact on the study of biological systems by exploring various structure–function relationships.[59] In the context of this Account, MD has been utilized to study the role of substrate dynamics in reactivity and product distributions. This is a promising approach to design novel catalytic agents where specific products are desired in an otherwise promiscuous binding pocket like that of NphB. The prohibitive cost of extensive sampling, however, poses a major barrier to utilizing this promising technique for large-scale studies aimed at this goal. QM/MM methods[32] allowed us to garner unique insights into enzymatic catalysis, and by the nature of this method we are typically focusing on a very specific reaction process where large-scale conformational changes are generally not encountered. Nonetheless, sampling especially using more accurate QM models will be an issue that will have to be further addressed in the coming years. Indeed, computational biology studies, at the molecular level, have two significant issues that will need to be further addressed in the future: We must accurately[60−65] calculate the energies and forces involved in these systems with very expensive QM models, while simultaneously sampling all relevant states of a system. Sophisticated QM models can address the accuracy issue, but how to extensively sample biological systems at the QM/MM level of theory using more accurate QM representations will be an active area of future research.