Advances in Understanding Enzyme Function by Way of Computational Investigations Correlated with Experimental Observations

1. Introduction

Enzymes are biological catalysts that play a pivotal role in the essential chemical reactions taking place in living organisms and in technological applications. For example, they are needed for providing cellular energy, repairing tissues, organs and cells, stimulating brain functions, digesting food etc. Research interest in the mechanisms of enzyme-catalyzed reactions and their kinetics dates from the middle of the 18th century on, particularly in terms of the relation between enzyme structure and the determinants of the specificity and efficiency of catalysis (Fersht, 1985). Experimental and biomolecular modeling studies, based on the lock-and-key, induced-fit and transition-state stabilization hypotheses, enabled to conceive enzymes as structural scaffolds that selectively recognize and bind substrates and cofactors (Fischer, 1894; Koshland, 1958; Robertus et al., 1972) in a reaction environment that is different from the bulk solvent (Bruice & Benkovic, 2000). They are able to speed up reactions by about 10¹⁷ folds or more (Richard, 2013). Active site residues that are in direct contact with substrates and cofactors are critical for orienting other reactants and keeping them in fixed positions properly (Nagel & Klinman, 2009; Nagel et al., 2013), indicating that mutations of these specific residues allow for the change of the rate of reactions in a desired way (Nagel et al., 2013). Active sites presumably have shapes that are complementary to the transition states of catalyzed reactions, and thus provide complementary structural and electrostatic environments for the preferential binding and/or stabilization of the transition states (Warshel et al., 2006). Electrostatic effects are one of the key factors for understanding the mechanism of enzyme catalysis. In contrast to the classical view that considers enzymes as rigid receptors responsible for interacting with properly-oriented ligands (Figure 1, left), growing evidences suggest essential impact of the bulk and hydration-shell solvent (Figure 1, right) on enzymatic structure-function relationships (Agarwal, 2005; Fenimore et al., 2004). The question of how other factors than direct structural interactions, such as long-range effects (electrostatic effects, external and internal solvent molecules) and the energy released from the binding of substrate and/or cofactor, contribute to catalysis is a research challenge. Noteworthy is that bio-organic and protein molecules, designed to mimic enzyme active site, are mainly associated with lower catalytic efficiencies than their natural counterparts (Benkovic, 1992; Siegel et al., 2010). In other words, a particular mechanism by which each factor contributes to overcoming the activation energy barrier depends on enzyme system and type of chemical catalysis.

Understanding enzymes as intrinsically flexible macromolecules implicated substantial interest in studying the impact of internal motions, driven by environment and temperature, on their function (Antoniou & Schwartz, 2011; Arora & Brooks, 2009; Aviram et al., 2018; Boehr et al., 2018; Choy et al., 2017; Clouthier et al., 2012; Garcia-Viloca et al., 2004; Henzler-Wildman et al., 2007; Hirschi et al., 2010; Holliday et al., 2017; Kim et al., 2017; Klinman, 2013; Klinman et al., 2017; Kohen, 2015; Kondo et al., 2017; Lisi & Loria, 2016; Machleder et al., 2010; Micheletti, 2013; Onions et al., 2016; Otten et al., 2018; Quaytman & Schwartz, 2007; Schramm et al., 2008). The role of internal motions is a matter of intensive debate in the literature with several contradictory views (Benkovic et al., 2008; Boehr et al., 2006; Kamerlin & Warshel, 2010; Kohen, 2015; Nagel & Klinman, 2009; Ramanathan & Agarwal, 2011) due to various definitions of protein motions (Kohen, 2015). It means that some particular conformational fluctuations may be associated with encoding dynamics and the catalytic rate acceleration being enzyme’s features that promote function (Agarwal, 2019). The origin and hierarchy of these internal motions in enzyme structure are discussed hereafter. Also, recent applications of hybrid quantum mechanics/molecular mechanics (QM/MM) calculations in understanding enzyme function, which are well calibrated relative to experimental observations by allowing the change of electronic structure in the functionally critical regions of enzyme structure, are reviewed.

2. Reaction-promoting conformational fluctuations in the enzyme structure

The term enzyme dynamics often has a misleading meaning, considering a wide range of motions associated with. A physically meaningful distinction of the internal motions is based upon their timescales (Agarwal, 2019). The first type of motions includes those at the sub-picosecond timescales, involving a small number of atoms or a few residues, sampling the same area of conformational landscape, averaging around the same mean structure over time, and being driven by temperature (Figure 2a). The second type of motions includes those at the nanosecond, microsecond and higher timescales, involving fluctuations within protein domains (loops, secondary structure elements) or entire protein, sampling remote areas of conformational landscape, and being driven by coupling global and local vibrational modes in dynamic allostery of proteins (Figure 2b). The third type of motions includes temporary (or even permanent) change in protein structure, shifting the thermodynamic equilibrium through changes in the energetic landscape due to the binding (or release) of other molecules (Figure 2c), and leading to a change in average structure (Agarwal, 2019; Hawkins & McLeish, 2006).

A variety of experimental methods, including nuclear magnetic resonance (NMR) (Boehr et al., 2006; Holliday et al., 2017; Palmer, 2015), single molecule techniques (Zhang et al., 2004), hydrogen/deuterium exchange (Zavodszky et al., 1998) and neutron scattering (Tehei et al., 2007), have been combined with computational methods (Garcia-Viloca et al., 2004; Holliday et al., 2017) in order to identify motions that are coupled to enzyme reaction mechanisms (Boehr et al., 2006; Klinman et al., 2017). It was shown that a number of enzymes, including cyclophilin A (Agarwal et al., 2004; Camilloni et al., 2014; Howard et al., 2003), dihydrofolate reductase (DHFR) (Agarwal et al., 2002a; Boehr et al., 2013; Sawaya & Kraut, 1997; Wang et al., 2006), liver alcohol dehydrogenase (Agarwal et al., 2000; Bandaria et al., 2008; Billeter et al., 2001; Nagel et al., 2013), ribonuclease A (Veeraraghavan et al., 2010), lipase B (Layfield & Hammes-Schiffer, 2013) and so on, contain functionally relevant residues being distant from the active site and associated with a hierarchy of conformational fluctuations that mediate the thermodynamic coupling of the surrounding solvent with enzyme catalysis (Agarwal, 2019). These particular residues are arranged in networks (Figure 3) linking loop regions on the enzyme surface to the active site (Agarwal, 2004; Agarwal et al., 2002b; Agarwal et al., 2004; Agarwal et al., 2012; Ramanathan & Agarwal, 2011) through internal motions having the rates that are measurable and reproducible in both apo enzyme and enzyme complexes (Boehr et al., 2006; Eisenmesser et al., 2002; Fraser et al., 2009). Thus, the nature of the protein flexibility is not random, but intrinsic to the topology of enzyme structure (Agarwal, 2019). Functionally important conformational substates can be leveraged to overcome the activation barrier through selectively directed energy transfer into the active site, while dynamic regions of the enzyme structure involved in slower conformational shifts are thought to be evolutionary conserved (Agarwal, 2019; Ramanathan & Agarwal, 2011). This biophysical standpoint stimulates the development of hyper-catalytic enzymes by identifying and modulating reaction-promoting conformational fluctuations (Agarwal et al., 2012).

Although allosteric modulation is a common mechanism for regulation of biochemical processes (Monod et al., 1963), an insufficiently clarified aspect of the control of enzyme activity in the presence of long-range effects was related to electrostatics and correlated structural interactions (Baldwin & Chothia, 1979) in terms of energy flow in the enzyme structure (Frauenfelder et al., 2001; Lee et al., 2008; Leitner, 2008; Reynolds et al., 2011). In addition to enzyme involvement in lowering the effective transition state barrier by raising the ground state and/or lowering the transition state through stabilization, additional barriers that need to be overcome by enzyme-catalyzed reactions were energetically estimated to cost circa 5-10 kcal mol^-1 (Agarwal, 2019). The random thermal fluctuations in enzyme structure provide about 1 kcal mol^-1 at room temperature. Thus, the remaining amount of energy of 4-9 kcal mol^-1 was hypothesized to come from the thermodynamic energy of the solvent on the protein surface (Agarwal, 2005; Agarwal, 2006; Ramanathan & Agarwal, 2011). It indicates that the pathways of connected residues, which span from the surface regions to the active site (Figure 3), serve as energetic channels for coupling the surrounding solvent to the enzyme reaction through intrinsic long-range effects (Fenimore et al., 2004). The thermodynamic coupling with solvent is a consequence of the energy produced by the fast motions of the surface loop residues with long side chains providing their high flexibility (Bellissent-Funel et al., 2016; Borreguero et al., 2011; Klinman et al., 2017; Ramanathan & Agarwal, 2009). This energy is transferred to the intermediate motions and eventually to the slow protein motions in the active site through preserved hydrogen bonds (Agarwal, 2005). This energy flow is not dissipative due to random collisions (Leitner, 2008; Moritsugu & Kidera, 2004), but likely directed flow of coupled conformational motions, mainly transfer of vibrational energy from a specific normal mode to several other resonant modes (Moritsugu & Kidera, 2004), as observed between α-helices (Leitner, 2008; Yu & Leitner, 2005). This aspect may be quantitatively rationalized using the configuartional entropy, being the entropy of the solute without translational and rotational contributions (Andricioaei & Karplus, 2001; Schlitter, 1993). The separation of translational entropy of the solute is generally acceptable without any approximation, while the removal of rotational entropy is associated with the assumption of negligible correlation between the internal degrees of freedom and the global rotation of the solute. The configurational entropy, estimated by the quasiharmonic approximation that treats a system of atoms with many degrees of freedom as a system of uncoupled harmonic oscillators and assumes a multivariate normal distribution of atomic fluctuations, gives the upper bound of the true entropy contribution to the free energy of non-covalent binding (Andricioaei & Karplus, 2001). The configurational entropy change can pragmatically be considered as the thermodynamic consequence of the change in receptor flexibility when binding ligands (Mitrasinovic, 2015, 2018, 2019). Understanding that slower global motions are consequences of faster local motions in proteins speaks in favor of the hypothesized model of enzyme activity (Figure 1, right) that was supported by a few experimental studies (Boehr et al., 2013; Chi et al., 2015; Eisenmesser et al., 2005; Gagne et al., 2012, 2015; Holliday et al., 2017; Wang et al., 2006).

Protein motions that occur at short timescales allow conformational sampling within the well (Figure 2a), while those at longer timescales stretch over multiple minima and allow interconversion between different conformational substates (Figure 2b). The rate of conformational transitions in proteins was related to slower anharmonic energy relaxation rate at the nanosecond timescale with the assistance of multilevel vibrational excitations (Xie et al., 2002). It was reported that functionally relevant conformational substates, characterized by small populations, contain structural and dynamic features that promote catalytic reactions (Ramanathan et al., 2014). The observed anharmonic nature of protein motions at long timescales (Borreguero et al., 2011; Ichiye & Karplus, 1988; Pontiggia et al., 2007) led to the development of a new computational method, termed “quasi anharmonic analysis” (QAA), for the identification of functionally relevant conformational substates by use of higher order statistics (Ramanathan et al., 2011). Conformational transitions obtained by QAA (Ramanathan et al., 2011) were in agreement with results of the NMR spin relaxation studies (Eisenmesser et al., 2002; Lange et al., 2008), showing that long-range fluctuations primarily refer to both loop areas on protein surface and loops located in the vicinity of the active site and represent rare, short-lived enzyme conformations that facilitate the acquisition of the transition states through their geometrical features (Agarwal, 2019). The impact of the hidden conformational substates on enzyme catalysis was highlighted experimentally (Fraser et al., 2009).

Although enzymes undergo a range of motions, the lack of clear definition of protein dynamics has been making a number of computational and experimental studies unsuccessful in providing consistent explanations on which particular protein motions promote catalysis (Agarwal, 2019). The development of models, based on the QAA methodology and followed by experimental validation, is a promising approach to getting mechanistic insights into the identification of functionally relevant conformational substates. For example, the study of enzyme activity in miscible organic solvents has shown that not all types of protein dynamics alter the catalytic rate, but a particular subset of conformational transitions into functionally relevant substates (Duff et al., 2018).

3. On the relevant structural basis of enzymatic function

Essential conformational fluctuations that subtly drive receptor structure into functionally relevant substates are reflected through the formation or break up of intra- and intermolecular interactions in the functionally critical regions of enzyme structure on the nanosecond timescale (Pavlov & Mitrasinovic, 2010). The advent of computers has enabled the development and application of theoretical models that aid in decoding ligand/enzyme binding as a major topic in medicine, chemistry and biology. All the main computational methods have some pros and cons, so that any particular choice should be carefully made after considering the chemical nature of the system under investigation, the available structural and experimental data, and the ultimate goal of the research (Limongelli, 2020).

The possibility of observing various conformations in static crystal structures, as a consequence of different crystallization conditions or procedures, implies that the dynamics of ligand/enzyme binding is a functionally remarkable standpoint of consideration. Theoretical studies that deal with the molecular dynamics (MD) of enzyme structure are based on molecular mechanics (MM) simulations associated with acceptable computational costs (Limongelli, 2020). This advantage is, however, balanced with the fact that enzyme electronic structure is treated as being unchanged throughout the entire simulation. Thus, it is hard to grasp the underlying charge transfer of intra- and intermolecular interactions by applying only MD methods (Pavlov & Mitrasinovic, 2010). Instead, it is necessary to exploit quantum mechanical (QM) methods that model the distributions of electrons in molecules explicitly. The geometries, electronic structure and reactions of small molecules are commonly studied by this approach. As an enzyme system typically consists of thousands of atoms means that the modeling of an entire enzyme structure using QM methods is a computationally challenging task. Furthermore, the need to account for solvent effects using adequate descriptions of the non-bonded interactions between solvent molecules, and between the solvent and solute, additionally makes required quantum chemical wave functions a lot more sophisticated. Therefore, fully realistic models of the solvation of macromolecular complexes are currently out of reach by calculation. In need of the strategies capable of treating efficiently charge-related effects in the functionally relevant regions of the large biomolecular structures, hybrid - QM/MM approach was developed with the ultimate goal to yield subtle conformational insights into the recognition of biological macromolecules. The general idea was to approximately treat electron correlation by QM in the key functional region of biomolecular structure, while having the rest of the system treated by MM (Lonsdale et al., 2012).

QM methods mutually differ depending on the approximations made to the Schrodinger equation that cannot be solved exactly for molecular systems containing more than one electron. Three main types of QM methods are semi-empirical approaches, ab initio and density functional theory (DFT). Semi-empirical molecular orbital methods, such as AM1 and PM3, are based on the Modified Neglect of Diatomic Differential Overlap (MNDO) approximation, the least computationally demanding and the least accurate of the QM methods. In some cases a MM treatment may be better than a low level semi-empirical one. Semi-empirical methods can be applied to larger systems than ab initio or DFT methods, typically to those of hundreds of atoms. The foundation of the most straightforward ab initio method is based on Hartree-Fock (HF) theory, which makes the very fundamental approximation that each electron's spatial distribution is independent of the other electrons’ instantaneous motion. By ignoring electron correlation or the tendency of electrons to avoid each other, HF calculations on chemical reactions are often misleading and impractical. Other ab initio methods that include electron correlation, such as those based on Moller-Plesset perturbation theory (e.g., MP2), configuration interaction (CI), or coupled cluster theory (CC), exploit HF wavefunctions as a starting point and provide a significant improvement in accuracy over HF calculations (Figure 4), but with much higher computational costs that make their application to systems with many tens of atoms difficult (Lonsdale et al., 2012). DFT provides a theoretical framework for determining the ground-state energy of a molecule from knowledge of the electron density distribution. The electron density is a function of only three variables, and thus much less complex for calculation than the ab intio wave function being a function of 3N variables, where N denotes the number of electrons. As an exact form of the functional for connecting the density to the energy is unknown, many approximate functionals have been empirically developed (Ahmed et al., 2012; Mitrasinovic, 2001, 2002, 2003a, 2003b, 2004, 2005, 2010a, 2012, 2013, 2014a, 2014b); e.g., B3LYP – a combination of the Becke’s three-parameter hybrid exchange functional (Becke, 1993) and the Lee-Yang-Parr correlation functional (Lee et al., 1988). DFT methods are capable of approaching accuracy of the correlated ab initio methods at much more reasonable computational expense (Lonsdale et al., 2012; Pavlov & Mitrasinovic, 2010), but often incapable of describing dispersive forces adequately (Grimme et al., 2007), particularly within unconventional systems (Hobza et al., 1995; Kristyan & Pulay, 1994). By including a dispersion correction (D) into DFT, DFT-D has become capable of performing the energy minimizations and vibrational analyses of extended molecular systems containing hundreds of atoms (Hachim et al., 2020; Jurecka et al., 2007; Mohammadkhani et al., 2020; Wang & Zhu, 2020). The second order Taylor expansion of the DFT total energy around a reference density has been a starting point for further approximations leading to the the Self-Consistent-Charge Density-Functional Tight-Binding (SCC-DFTB) model that contains additional induced-charges term for a proper description of charge-transfer phenomena (Elstner et al., 1998, 2001). SCC-DFTB with an empirical correction for dispersion energy (SCC-DFTB-D) (Grimme et al., 2007) has contributed to a more balanced description of the interactions inside the systems (Kubar et al., 2007). By being several orders of magnitude faster than DFT-D (Pavlov & Mitrasinovic, 2010) and providing geometries and relative energies comparable to DFT and ab initio calculations (Elstner et al., 2003), SCC-DFTB-D is an optimal, quantum chemical choice for modeling the structural parts of enzyme systems that are mainly responsible for their functions (Figure 4).

In order to provide insights into enzyme function by biomolecular modeling, a high-resolution (less than 2 Å) X-ray crystallographic structure of an enzyme complex is usually a starting point that tends to accurately represent the reacting system. A crystal structure is a result of averaging over all the conformations in the crystal and over the entire time of data acquisition. It means that the alternative conformations may be observed for some structural elements of the system, particularly for those that are very mobile and have no well-defined single conformation or position such as amino acid side chains, surface loops, and terminal regions. It is difficult to structurally isolate the active enzyme-substrate complexes as their reactions proceed too fast. The crystal structures of less efficient mutants or substrates can sometimes be resolved by varying redox conditions, if the reactions are slow enough to allow the observation of the complexes. It is widely accepted that the crystallographic structure of an enzyme-inhibitor complex, recognized as a closed or “receptor inaccessible” conformation, is an appropriate choice, where the inhibitor mimics the substrate, product, and transition state (or an intermediate) in its bound conformation (Lonsdale et al., 2012). Therefore, an apo (ligand-free) crystal structure of an enzyme, denoted by an open or “receptor accessible” conformation, may be seen as one of little use, by being incapable of predicting conformational changes upon binding. However, it has been shown that the open conformations in the structures of several enzymes, such as influenza A/H5N1 virus neuraminidase (Mihajlovic & Mitrasinovic, 2008, 2009a, 2009b; Mitrasinovic, 2009, 2010b; Wang et al., 2010), HIV-1 protease (Yu et al., 2017), Yersinia protein tyrosine phosphatase (YopH) (Mitrasinovic, 2023) and others, should be considered as conformations of some relevance for exploring the dynamics of their recognition by ligands. In all those cases, a common point is that an open state is defined by the conformation of a loop that is located nearby the active site and affects ligand binding through the conformational flexibility of some of its residues (Figure 2c). For example, the apo HIV-1 protease shows a variety of flap conformations, ranging from closed (Pillai et al., 2001; Robbins et al., 2010) to semi open (Heaslet et al., 2007; Lapatto et al., 1989; Liu et al., 2006; Spinelli et al., 1991) and to wide open (Coman et al., 2008; Heaslet et al., 2007; Martin et al., 2005). It means that the wild-type (wt) receptor may initially approach to the ligand in an open conformation and, then, likely adopts the higher energy state or closed conformation of the particular loop via a subtle conformational change. A more realistic interpretation, based on the thermodynamics of protein-ligand binding, is that a semi-open intermediate reduces the free energy barrier between the closed and open states (Gur et al., 2020), despite the fact that it is difficult to capture a semi-open state by currently available experimental techniques (Lapatto et al., 1989; Spinelli et al., 1991; Navia et al., 1989) due to its short transient lifetime. In the case of YopH, this intermediate has been observed by means of computational simulations, which demonstrate that the closure of the so-called “WPD” loop occurs in a half-closed state in the apo, native protein structure (Hu & Erec Stebbins, 2006). The possibility of having the conformational flexibility of the WPD loop correlated with the conformational behavior of another loop as the second binding site in the YopH structure (Hu & Erec Stebbins, 2006) has been elaborated in more detail (Mitrasinovic, 2024a, 2024b), thereby indicating that the potential allosteric binding sites offer new opportunities for the design of novel, selective YopH inhibitors (Mitrasinovic, 2023, 2024a, 2024b) as antibacterial agents (Kuban-Jankowska et al., 2016). It has also been reported that the catalytic process during glycolysis is associated with the synergistic effect of conformational conversions through a semi-open state in phosphoglycerate kinase 1 (PGK1) due to the releasing/binding of products/substrates, as well as with post-translational modifications, thereby enabling to better understand the PGK1 substrates/products conversion process and the development of small drug molecules targeting PGK1 (Liu et al., 2022). A similar approach would apply to the structural understanding of the inhibition mechanism of a lipid peroxidase pathway associated with therapy-resistant cancer cells in various mesenchymal cell state contexts (Viswanathan et al., 2017). It has even been suggested that forming the semi-open or intermediate state in a protein structure is a prerequisite to the recognition of host cells (Gur et al., 2020). In some cases, specific single point mutations might affect the attachment of proteins to cell surface receptors by improving the stability of proteins in their open and closed conformations and driving the thermodynamic preferences toward the open mutant form (Verkhivker et al., 2022).

Cytochrome P450 enzymes (P450s or CYPs) are mainly involved in Phase I drug metabolism, which includes the introduction of polar functional groups into non-polar molecules and can lead to inactivation or activation of drugs. Therefore, modeling P450-catalyzed reactions is important for drug discovery (Lonsdale et al., 2012). QM calculations using small cluster models improved the understanding of the mechanisms through which these enzymes stimulate oxidative metabolism (Bathelt et al., 2003, 2004, 2008; Guallar et al., 2003; Shaik et al., 2005, 2010; Zurek et al., 2006). The modeling of such systems, based on the B3LYP density functional, is in good agreement with the available experimental data (Bathelt et al., 2008; Harvey et al., 2006; Lonsdale et al., 2010, 2011; Shaik et al., 2010; Zurek et al., 2006). The determination of the crystal structures of P450 enzymes, combined with the development of QM (DFT)/MM methods (Figure 4), made it possible to model the mechanisms of oxidation in the enzymatic environment. The QM/MM (SCC-DFTB-D/Amber-FF99SB) methodology, being well calibrated against relevant experimental data, was an inevitable choice between alternatives to address several vital issues by studying the ligand/enzyme binding interactions and evaluating the competitive inhibition parameters. Both quantitative and qualitative insights into the functional ability of structurally diversified flavonoids to inhibit P4503A4, the most versatile enzyme involved in drug metabolism, were reported (Mitrasinovic, 2022). The same approach was employed to rationalize the affinity and specificity of peptide inhibitors to caspase-2, one of the most evolutionarily conserved caspases that is involved in cell-cycle regulation and DNA repair as well as in mediating neurodegenerative disorders and non-alcoholic steatohepatitis pathogenesis (Mitrasinovic, 2020). In a similar manner, the active/inactive functional behavior of YopH, an extraordinarily active protein tyrosine phosphatase involved in breaking down signal transduction mechanisms in immune cells and inhibiting the immune response, was investigated by observing its modes of recognition by carboxylic acid derivatives (Mitrasinovic, 2023) and identifying a pharmacophore for designing novel inhibitors capable of achieving specificity against phosphatases (Mitrasinovic, 2024a, 2024b). In all the QM/MM models (Mitrasinovic, 2020, 2022, 2023), the QM region was defined to comprise not only active site residues and the ligand, but also loops proximal to the active site and flexible surface loop regions previously shown by using experimental means to promote function through their conformational fluctuations within the enzymatic structure of P4503A4 (Sevrioukova & Poulos, 2015), caspase-2 (Tang et al., 2011), and YopH (Schubert et al., 1995), respectively. Thus, SCC-DFTB-D essentially reconciled an implicit description of charge-transfer phenomena with a balanced description of the interactions inside the systems, enabling to quantitatively express the impact of the underlying conformational substates on the enzyme’s performances. It was for this reason necessary to avoid unrealistic effects by choosing the adequate positions of atoms on the QM/MM frontier (Figure 4), so that they do not cut across any polar bond (Mitrasinovic, 2020, 2022, 2023). A pseudo-QM atom can be added to create an appropriate bond, e.g., modeling edge atoms as hydrogen atoms and neglecting some of their interactions with MM atoms. The charges of MM groups next to QM groups are commonly set to zero in order to escape unbalanced charge interactions across the boundary, even though more sophisticated electrostatic settings can be applied too (Lonsdale et al., 2012). Finally, in order to reasonably describe conformational fluctuations within an enzyme structure as pathways of energy transfer that encode dynamics, an evaluation of chosen protein force field for MM part is needed in terms of how accurately it predicts the structure of multiple flexible loop domains present in proteins (Love et al., 2022).

To sum up, it is pragmatic to note two sides - the mechanistic and electronic structure aspects of modeling enzyme-catalyzed reactions. The mechanistic aspect of this topic can be explored using the QM/MM strategy, the empirical valence bond (EVB) methodology and the cluster approach. Although useful, the cluster approach is out of the scope of this review, because it includes only a small portion of the enzyme structure in the models for exploring chemical properties and reaction mechanisms (Himo, 2017). An insufficient amount of available experimental data about the catalytic efficiencies of enzyme substrates often means that a more complete understanding of the true substrate preferences is hardly possible. The above-mentioned QM/MM results (Mitrasinovic, 2020, 2022, 2023) illustrate the way in which additional experiments can be stimulated by biomolecular modeling, in order to better understand the kinetics of ligand/enzyme recognition (e.g., in terms of specificity (Mitrasinovic, 2024a, 2024b)) for further research or therapeutic product development (Mitrasinovic, 2020). We advocate that both mechanistic and electronic structure standpoints can be reconciled by applying a more advanced strategy based on molecular dynamics simulations with the QM/MM potentials (Figure 4). For example, in the case of hydrolysis by glutamate carboxypeptidase, it was possible to construct reaction energy profiles with realistic barrier heights on the reaction pathways for the wt and mutated enzymes (Krivitskaya et al., 2021). The QM/MM potentials were evaluated by means of atom-centered basis sets and various hybrid DFT functionals in QM parts, as well as commonly accepted force field parameters in MM parts (Krivitskaya et al., 2021). To generate MD trajectories on 10 picosecond timescale (Krivitskaya et al., 2021), the umbrella sampling (US) scan of the Gibbs free energy landscape was exploited with the subsequent statistical analysis (Limongelli, 2020), using the umbrella integration methods to quantify the elementary steps of chemical reactions in the enzyme active sites (Krivitskaya et al., 2021). Anyway, the conformational flexibility of molecules (large proteins or small ligands) is an important feature. It means that conformational sampling that generates an ensemble of biologically meaningful conformations is a precondition for computer-aided drug design. It was observed that funnel metadynamics (FM) at the nanosecond timescales leads to an efficient characterization of the binding free energy surface and an accurate calculation of the absolute ligand/protein (or enzyme) binding free energy compared to an experimental value (Limongelli et al., 2013). FM also yields important insights into the binding process, such as the presence of alternative binding modes and the role of water, at an affordable computational cost (Limongelli, 2020; Limongelli et al., 2013). Some of the generally applicable correction terms can further improve the accuracy of physical pathways methods for calculating the free energy, such as Markov state modeling, FM, US and steered MD (Hall et al., 2020). Continued advances in increasing computational power are expected to make QM/MM MD calculations (Figure 4) routine in the timescales involving conformational fluctuations within protein domains.

4. Future outlook

The origin of the enormous catalytic power of enzymes is one of the key questions in biology today. To address this issue, besides electrostatic effects that stabilize the transition states of enzymatic reactions, it is necessary to investigate the impact of other possible factors (Pisliakov et al., 2009). An interesting possibility is that the dynamics of conformational transitions makes a remarkable contribution to enzyme catalysis through a hierarchy of internal motions in intrinsically flexible enzyme structures (Agarwal, 2019). A multiscale approach has been developed to explore the dynamic nature of enzyme catalysis in the millisecond timescale, indicating that there exists basically no coupling between chemistry and protein structural fluctuations as long as the chemical barrier is greater than a few kcal/mol (which is the case in enzymatic reactions) (Pisliakov et al., 2009). In spite of that, the nature of this coupling is a question of great importance that deserves continued research interest. This view may be relevant for the development of hyper-catalytic enzymes by identifying and modulating reaction-favoring conformational fluctuations (Agarwal, 2019; Agarwal et al., 2012).

Advanced QM/MM MD calculations are getting feasible due to advances in computer power that enables to reconcile mechanistic and electronic structure aspects (two sides of the same enzyme-catalyzed reaction) in the nanosecond, microsecond and higher timescales (Figure 4). In the context of ligand/protein (or enzyme) recognition, advanced MD techniques, such as replica exchange methods (Yang et al., 2015; Yang & MacKerell, 2015), can enhance the sampling efficiency of the conformational landscape at the microsecond scales and reveal new modes of interaction (Mitrasinovic, 2019). The integration of more rigorous physics-based methods, such as emerging methods for ensemble based virtual screening, enables the use of multiple enzyme conformations, reveals new modes of activation/inhibition, provides a better understanding of the structural dynamics of disease-related target enzymes, and supports the development of novel patient-specific therapeutics that are less prone to drug resistance (Amaro & Li, 2010). Thus, computational protocols, such as SILCS (site identification by ligand competitive saturation) for structure-based drug design and CSP (conformationally sampled pharmacophore) for ligand-based drug design, are more advantageous in practice than CADD (computer-aided drug design) methods being based on single crystal structure or a limited number of ligand conformations (Yu & MacKerell, 2017).