Exploring the Minimum-Energy Pathways and Free-Energy Profiles of Enzymatic Reactions with QM/MM Calculations

High-performance QM/MM Enhanced Sampling Molecular Dynamics Simulations with GENESIS SPDYN and QSimulate-QM

A new module for quantum mechanical/molecular mechanical (QM/MM) calculations is implemented in a molecular dynamics (MD) program, SPDYN in GENESIS, interfaced with an electronic structure program, QSimulate-QM. The periodic boundary condition (PBC) in QM/MM simulation is incorporated via QM calculation in real space with duplicated MM charges and particle mesh Ewald (PME) calculation with QM and MM charges. A highly optimized code is implemented in QSimulate-QM, particularly for the density functional tight-binding (DFTB) method, where the interaction between the QM and MM regions is computed utilizing multipole expansions. Together with highly parallelized algorithms in SPDYN, the developed program performs MD simulations based on DFTB in the QM size of ∼100 atoms and MM of ∼100,000 atoms with a better performance than 1 ns/day using one computer node. This feature paves the way for QM/MM-MD enhanced sampling simulations. Various enhanced sampling methods in GENESIS, namely, generalized replica exchange solute tempering (gREST), replica-exchange umbrella sampling (REUS), and path sampling with the string method, are demonstrated at the QM/MM level to compute the Ramachandran plot of alanine dipeptide and the potential of mean force (PMF) of a proton transfer reaction in an enzyme.

A structural perspective on enzymes and their catalytic mechanisms

In this perspective, we analyse the progress made in our knowledge of enzyme sequences, structures and functions in the last 2 years. We review how much new enzyme data have been garnered and annotated, derived from the study of proteins using structural and computational approaches. Recent advances towards capturing 'Catalysis in silico' are described, including knowledge and predictions of enzyme structures, their interactions and mechanisms. We highlight the flood of enzyme data, driven by metagenomic sequencing, the improved enzyme data resources, the high coverage in Protein Data Bank of E.C. classes and the AI-driven structure prediction techniques that facilitate the accurate prediction of protein structures. We note the focus on disordered regions in the context of enzyme regulation and specificity and comment on emerging bioinformatic approaches that capture reaction mechanisms computationally for comparing and predicting enzyme mechanisms. We also consider the drivers of progress in this field in the next five years.

Efficient Training of Neural Network Potentials for Chemical and Enzymatic Reactions by Continual Learning

Machine learning (ML) methods have emerged as an efficient surrogate for high-level electronic structure theory, offering precision and computational efficiency. However, the vast conformational and chemical space remains challenging when constructing a general force field. Training data sets typically cover only a limited region of this space, resulting in poor extrapolation performance. Traditional strategies must address this problem by training models from scratch using old and new data sets. In addition, model transferability is crucial for general force field construction. Existing ML force fields, designed for closed systems with no external environmental potential, exhibit limited transferability to complex condensed phase systems such as enzymatic reactions, resulting in inferior performance and high memory costs. Our ML/MM model, based on the Taylor expansion of the electrostatic operator, showed high transferability between reactions in several simple solvents. This work extends the strategy to enzymatic reactions to explore the transferability between more complex heterogeneous environments. In addition, we also apply continual learning strategies based on memory data sets to enable autonomous and on-the-fly training on a continuous stream of new data. By combining these two methods, we can efficiently construct a force field that can be applied to chemical reactions in various environmental media.

Multiscale calculations reveal new insights into the reaction mechanism between KRASG12C and α, β-unsaturated carbonyl of covalent inhibitors

Utilizing α,β-unsaturated carbonyl group as Michael acceptors to react with thiols represents a successful strategy for developing KRASG12C inhibitors. Despite this, the precise reaction mechanism between KRASG12C and covalent inhibitors remains a subject of debate, primarily due to the absence of an appropriate residue capable of deprotonating the cysteine thiol as a base. To uncover this reaction mechanism, we first discussed the chemical reaction mechanism in solvent conditions via density functional theory (DFT) calculation. Based on this, we then proposed and validated the enzymatic reaction mechanism by employing quantum mechanics/molecular mechanics (QM/MM) calculation. Our QM/MM analysis suggests that, in biological conditions, proton transfer and nucleophilic addition may proceed through a concerted process to form an enolate intermediate, bypassing the need for a base catalyst. This proposed mechanism differs from previous findings. Following the formation of the enolate intermediate, solvent-assisted tautomerization results in the final product. Our calculations indicate that solvent-assisted tautomerization is the rate-limiting step in the catalytic cycle under biological conditions. On the basis of this reaction mechanism, the calculated kinact/ki for two inhibitors is consistent well with the experimental results. Our findings provide new insights into the reaction mechanism between the cysteine of KRASG12C and the covalent inhibitors and may provide valuable information for designing effective covalent inhibitors targeting KRASG12C and other similar targets.

Condensed Matter Systems Exposed to Radiation: Multiscale Theory, Simulations, and Experiment

This roadmap reviews the new, highly interdisciplinary research field studying the behavior of condensed matter systems exposed to radiation. The Review highlights several recent advances in the field and provides a roadmap for the development of the field over the next decade. Condensed matter systems exposed to radiation can be inorganic, organic, or biological, finite or infinite, composed of different molecular species or materials, exist in different phases, and operate under different thermodynamic conditions. Many of the key phenomena related to the behavior of irradiated systems are very similar and can be understood based on the same fundamental theoretical principles and computational approaches. The multiscale nature of such phenomena requires the quantitative description of the radiation-induced effects occurring at different spatial and temporal scales, ranging from the atomic to the macroscopic, and the interlinks between such descriptions. The multiscale nature of the effects and the similarity of their manifestation in systems of different origins necessarily bring together different disciplines, such as physics, chemistry, biology, materials science, nanoscience, and biomedical research, demonstrating the numerous interlinks and commonalities between them. This research field is highly relevant to many novel and emerging technologies and medical applications.

QM/MM Free Energy Calculations of Long-Range Biological Protonation Dynamics by Adaptive and Focused Sampling

Water-mediated proton transfer reactions are central for catalytic processes in a wide range of biochemical systems, ranging from biological energy conversion to chemical transformations in the metabolism. Yet, the accurate computational treatment of such complex biochemical reactions is highly challenging and requires the application of multiscale methods, in particular hybrid quantum/classical (QM/MM) approaches combined with free energy simulations. Here, we combine the unique exploration power of new advanced sampling methods with density functional theory (DFT)-based QM/MM free energy methods for multiscale simulations of long-range protonation dynamics in biological systems. In this regard, we show that combining multiple walkers/well-tempered metadynamics with an extended system adaptive biasing force method (MWE) provides a powerful approach for exploration of water-mediated proton transfer reactions in complex biochemical systems. We compare and combine the MWE method also with QM/MM umbrella sampling and explore the sampling of the free energy landscape with both geometric (linear combination of proton transfer distances) and physical (center of excess charge) reaction coordinates and show how these affect the convergence of the potential of mean force (PMF) and the activation free energy. We find that the QM/MM-MWE method can efficiently explore both direct and water-mediated proton transfer pathways together with forward and reverse hole transfer mechanisms in the highly complex proton channel of respiratory Complex I, while the QM/MM-US approach shows a systematic convergence of selected long-range proton transfer pathways. In this regard, we show that the PMF along multiple proton transfer pathways is recovered by combining the strengths of both approaches in a QM/MM-MWE/focused US (FUS) scheme and reveals new mechanistic insight into the proton transfer principles of Complex I. Our findings provide a promising basis for the quantitative multiscale simulations of long-range proton transfer reactions in biological systems.

Computational Insights into Amide Bond Formation Catalyzed by the Condensation Domain of Nonribosomal Peptide Synthetases

Nonribosomal peptide synthetases (NRPSs) are important enzymes that synthesize an array of nongenetically encoded peptides. The latter have diverse physicochemical properties and roles. NRPSs are modular enzymes in which, for example, the condensation (C-) domain catalyzes the formation of amide bonds. The NRPS tyrocidine synthetase from Brevibacillus brevis is responsible for synthesizing the cyclic-peptide antibiotic tyrocidine. The first step is formation of an amide bond between a proline and phenylalanine which is catalyzed by a C-domain. In this study, a multiscale computational approach (molecular dynamics and QM/MM) has been used to investigate substrate binding and catalytic mechanism of the C-domain of tyrocidine synthetase. Overall, the mechanism is found to proceed through three exergonic steps in which an active site Histidine, His222, acts as a base and acid. First, His222 acts as a base to facilitate nucleophilic attack of the prolyl nitrogen at the phenylalanyl's carbonyl carbon. This is also the rate-limiting step with a free energy barrier of 38.8 kJ mol-1. The second step is collapse of the resulting tetrahedral intermediate with cleavage of the S-C bond between the phenylalanyl and its Ppant arm, along with formation of the above amide bond. Meanwhile, the now protonated His222 imidazole has rotated toward the newly formed thiolate of the Ppant arm. In the final step, His222 acts as an acid, protonating the thiolate and regenerating a neutral His222. The overall mechanism is found to be exergonic with the final product complex being 46.3 kJ mol-1 lower in energy than the initial reactant complex.

GENESIS 2.1: High-Performance Molecular Dynamics Software for Enhanced Sampling and Free-Energy Calculations for Atomistic, Coarse-Grained, and Quantum Mechanics/Molecular Mechanics Models

GENeralized-Ensemble SImulation System (GENESIS) is a molecular dynamics (MD) software developed to simulate the conformational dynamics of a single biomolecule, as well as molecular interactions in large biomolecular assemblies and between multiple biomolecules in cellular environments. To achieve the latter purpose, the earlier versions of GENESIS emphasized high performance in atomistic MD simulations on massively parallel supercomputers, with or without graphics processing units (GPUs). Here, we implemented multiscale MD simulations that include atomistic, coarse-grained, and hybrid quantum mechanics/molecular mechanics (QM/MM) calculations. They demonstrate high performance and are integrated with enhanced conformational sampling algorithms and free-energy calculations without using external programs except for the QM programs. In this article, we review new functions, molecular models, and other essential features in GENESIS version 2.1 and discuss ongoing developments for future releases.

Learning QM/MM potential using equivariant multiscale model

The machine learning (ML) method emerges as an efficient and precise surrogate model for high-level electronic structure theory. Its application has been limited to closed chemical systems without considering external potentials from the surrounding environment. To address this limitation and incorporate the influence of external potentials, polarization effects, and long-range interactions between a chemical system and its environment, the first two terms of the Taylor expansion of an electrostatic operator have been used as extra input to the existing ML model to represent the electrostatic environments. However, high-order electrostatic interaction is often essential to account for external potentials from the environment. The existing models based only on invariant features cannot capture significant distribution patterns of the external potentials. Here, we propose a novel ML model that includes high-order terms of the Taylor expansion of an electrostatic operator and uses an equivariant model, which can generate a high-order tensor covariant with rotations as a base model. Therefore, we can use the multipole-expansion equation to derive a useful representation by accounting for polarization and intermolecular interaction. Moreover, to deal with long-range interactions, we follow the same strategy adopted to derive long-range interactions between a target system and its environment media. Our model achieves higher prediction accuracy and transferability among various environment media with these modifications.

Perspectives on Computational Enzyme Modeling: From Mechanisms to Design and Drug Development

Understanding enzyme mechanisms is essential for unraveling the complex molecular machinery of life. In this review, we survey the field of computational enzymology, highlighting key principles governing enzyme mechanisms and discussing ongoing challenges and promising advances. Over the years, computer simulations have become indispensable in the study of enzyme mechanisms, with the integration of experimental and computational exploration now established as a holistic approach to gain deep insights into enzymatic catalysis. Numerous studies have demonstrated the power of computer simulations in characterizing reaction pathways, transition states, substrate selectivity, product distribution, and dynamic conformational changes for various enzymes. Nevertheless, significant challenges remain in investigating the mechanisms of complex multistep reactions, large-scale conformational changes, and allosteric regulation. Beyond mechanistic studies, computational enzyme modeling has emerged as an essential tool for computer-aided enzyme design and the rational discovery of covalent drugs for targeted therapies. Overall, enzyme design/engineering and covalent drug development can greatly benefit from our understanding of the detailed mechanisms of enzymes, such as protein dynamics, entropy contributions, and allostery, as revealed by computational studies. Such a convergence of different research approaches is expected to continue, creating synergies in enzyme research. This review, by outlining the ever-expanding field of enzyme research, aims to provide guidance for future research directions and facilitate new developments in this important and evolving field.

Unraveling the catalytic mechanism of SARS-CoV-2 papain-like protease with allosteric modulation of C270 mutation using multiscale computational approaches

Papain-like protease (PL^pro) is a promising therapeutic target against SARS-CoV-2, but its restricted S1/S2 subsites pose an obstacle in developing active site-directed inhibitors. We have recently identified C270 as a novel covalent allosteric site for SARS-CoV-2 PL^pro inhibitors. Here we present a theoretical investigation of the proteolysis reaction catalyzed by the wild-type SARS-CoV-2 PL^pro as well as the C270R mutant. Enhanced sampling MD simulations were first performed to explore the influence of C270R mutation on the protease dynamics, and sampled thermodynamically favorable conformations were then submitted to MM/PBSA and QM/MM MD simulations for thorough characterization of the protease-substrate binding and covalent reactions. The disclosed proteolysis mechanism of PL^pro, as characterized by the occurrence of proton transfer from the catalytic C111 to H272 prior to the substrate binding and with deacylation being the rate-determining step of the whole proteolysis process, is not completely identical to that of the 3C-like protease, another key cysteine protease of coronaviruses. The C270R mutation alters the structural dynamics of the BL2 loop that indirectly impairs the catalytic function of H272 and reduces the binding of the substrate with the protease, ultimately showing an inhibitory effect on PL^pro. Together, these results provide a comprehensive understanding at the atomic level of the key aspects of SARS-CoV-2 PL^pro proteolysis, including the catalytic activity allosterically regulated by C270 modification, which is crucial to the follow-up inhibitor design and development.

Hybrid Quantum Mechanical/Molecular Mechanical Methods For Studying Energy Transduction in Biomolecular Machines

Hybrid quantum mechanical/molecular mechanical (QM/MM) methods have become indispensable tools for the study of biomolecules. In this article, we briefly review the basic methodological details of QM/MM approaches and discuss their applications to various energy transduction problems in biomolecular machines, such as long-range proton transports, fast electron transfers, and mechanochemical coupling. We highlight the particular importance for these applications of balancing computational efficiency and accuracy. Using several recent examples, we illustrate the value and limitations of QM/MM methodologies for both ground and excited states, as well as strategies for calibrating them in specific applications. We conclude with brief comments on several areas that can benefit from further efforts to make QM/MM analyses more quantitative and applicable to increasingly complex biological problems.

Allosteric regulation of β-reaction stage I in tryptophan synthase upon the α-ligand binding

Tryptophan synthase (TRPS) is a bifunctional enzyme consisting of α- and β-subunits that catalyzes the last two steps of L-tryptophan (L-Trp) biosynthesis. The first stage of the reaction at the β-subunit is called β-reaction stage I, which converts the β-ligand from an internal aldimine [E(Ain)] to an α-aminoacrylate [E(A-A)] intermediate. The activity is known to increase 3-10-fold upon the binding of 3-indole-D-glycerol-3'-phosphate (IGP) at the α-subunit. The effect of α-ligand binding on β-reaction stage I at the distal β-active site is not well understood despite the abundant structural information available for TRPS. Here, we investigate the β-reaction stage I by carrying out minimum-energy pathway searches based on a hybrid quantum mechanics/molecular mechanics (QM/MM) model. The free-energy differences along the pathway are also examined using QM/MM umbrella sampling simulations with QM calculations at the B3LYP-D3/aug-cc-pVDZ level of theory. Our simulations suggest that the sidechain orientation of βD305 near the β-ligand likely plays an essential role in the allosteric regulation: a hydrogen bond is formed between βD305 and the β-ligand in the absence of the α-ligand, prohibiting a smooth rotation of the hydroxyl group in the quinonoid intermediate, whereas the dihedral angle rotates smoothly after the hydrogen bond is switched from βD305-β-ligand to βD305-βR141. This switch could occur upon the IGP-binding at the α-subunit, as evidenced by the existing TRPS crystal structures.

The Impact of DFT Functional, Cluster Model Size, and Implicit Solvation on the Structural Description of Single-Metal-Mediated DNA Phosphodiester Bond Cleavage: The Case Study of APE1

Phosphodiester bond hydrolysis in nucleic acids is a ubiquitous reaction that can be facilitated by enzymes called nucleases, which often use metal ions to achieve catalytic function. While a two-metal-mediated pathway has been well established for many enzymes, there is growing support that some enzymes require only one metal for the catalytic step. Using human apurinic/apyrimidinic endonuclease (APE1) as a prototypical example and cluster models, this study clarifies the impact of DFT functional, cluster model size, and implicit solvation on single-metal-mediated phosphodiester bond cleavage and provides insight into how to efficiently model this chemistry. Initially, a model containing 69 atoms built from a high-resolution X-ray crystal structure is used to explore the reaction pathway mapped by a range of DFT functionals and basis sets, which provides support for the use of standard functionals (M06-2X and B3LYP-D3) to study this reaction. Subsequently, systematically increasing the model size to 185 atoms by including additional amino acids and altering residue truncation points highlights that small models containing only a few amino acids or β carbon truncation points introduce model strains and lead to incorrect metal coordination. Indeed, a model that contains all key residues (general base and acid, residues that stabilize the substrate, and amino acids that maintain the metal coordination) is required for an accurate structural depiction of the one-metal-mediated phosphodiester bond hydrolysis by APE1, which results in 185 atoms. The additional inclusion of the broader enzyme environment through continuum solvation models has negligible effects. The insights gained in the present work can be used to direct future computational studies of other one-metal-dependent nucleases to provide a greater understanding of how nature achieves this difficult chemistry.

Use of multistate Bennett acceptance ratio method for free-energy calculations from enhanced sampling and free-energy perturbation

Multistate Bennett acceptance ratio (MBAR) works as a method to analyze molecular dynamics (MD) simulation data after the simulations have been finished. It is widely used to estimate free-energy changes between different states and averaged properties at the states of interest. MBAR allows us to treat a wide range of states from those at different temperature/pressure to those with different model parameters. Due to the broad applicability, the MBAR equations are rather difficult to apply for free-energy calculations using different types of MD simulations including enhanced conformational sampling methods and free-energy perturbation. In this review, we first summarize the basic theory of the MBAR equations and categorize the representative usages into the following four: (i) perturbation, (ii) scaling, (iii) accumulation, and (iv) full potential energy. For each, we explain how to prepare input data using MD simulation trajectories for solving the MBAR equations. MBAR is also useful to estimate reliable free-energy differences using MD trajectories based on a semi-empirical quantum mechanics/molecular mechanics (QM/MM) model and ab initio QM/MM energy calculations on the MD snapshots. We also explain how to use the MBAR software in the GENESIS package, which we call mbar_analysis, for the four representative cases. The proposed estimations of free-energy changes and thermodynamic averages are effective and useful for various biomolecular systems.

Finite temperature string method with umbrella sampling using path collective variables: application to secondary structure change in a protein

The transition of an α-helix to a β-sheet in proteins is among the most complex conformational changes seen in biomolecular systems. Due to long time scales involved in the transition, it is challenging to study such protein conformational changes using direct molecular dynamics simulations. This limitation is typically overcome using an indirect approach wherein one computes the free energy landscape associated with the transition. Computation of free energy landscapes, however, requires a suitable set of collective variables that describe the transition. In this work, we demonstrate the use of path collective variables [D. Branduardi, F. L. Gervasio and M. Parrinello, J. Chem. Phys., 2007, 126, 054103] and combine it with the finite temperature string (FTS) method [E. Weinan, W. Ren and E. Vanden-Eijnden, J. Phys. Chem. B, 2005, 109, 6688-6693] to determine the molecular mechanisms involved during the structural transition of the mini G-protein from an α-helix to a β-hairpin. The transition from the α-helix proceeds via unfolding of the terminal residues, giving rise to a β-turn unfolded intermediate to eventually form the β-hairpin. Our proposed algorithm uses umbrella sampling simulations to simulate images along the string and the weighted histogram analysis to compute the free energy along the computed transition path. This work demonstrates that the string method in combination with path collective variables can be exploited to study complex protein conformational changes such as a complete change in the secondary structure.

Computer-aided drug design, quantum-mechanical methods for biological problems

Quantum chemistry enables to study systems with chemical accuracy (<1 kcal/mol from experiment) but is restricted to a handful of atoms due to its computational expense. This has led to ongoing interest to optimize and simplify these methods while retaining accuracy. Implementing quantum mechanical (QM) methods on modern hardware such as multiple-GPUs is one example of how the field is optimizing performance. Multiscale approaches like the so-called QM/molecular mechanical method are gaining popularity in drug discovery because they focus the application of QM methods on the region of choice (e.g., the binding site), while using efficient MM models to represent less relevant areas. The creation of simplified QM methods is another example, including the use of machine learning to create ultra-fast and accurate QM models. Herein, we summarize recent advancements in the development of optimized QM methods that enhance our ability to use these methods in computer aided drug discovery.

Engineering and Design of Programmable Genome Editors

Programmable genome editors are enzymes that can be targeted to a specific location in the genome for making site-specific alterations or deletions. The engineering, design, and development of sequence-specific editors has resulted in a dramatic increase in the precision of editing for nucleotide sequences. These editors can target specific locations in a genome, in vivo. The genome editors are being deployed for the development of genetically modified organisms for agriculture and industry, and for gene therapy of inherited human genetic disorders, cancer, and immunotherapy. Experimental and computational studies of structure, binding, activity, dynamics, and folding, reviewed here, have provided valuable insights that have the potential for increasing the functional efficiency of these gene/genome editors. Biochemical and biophysical studies of the specificities of natural and engineered genome editors reveal that increased binding affinity can be detrimental because of the increase of off-target effects and that the engineering and design of genome editors with higher specificity may require modulation and control of the conformational dynamics.

Mechanistic Insights into Enzyme Catalysis from Explaining Machine-Learned Quantum Mechanical and Molecular Mechanical Minimum Energy Pathways

With the increasing popularity of machine learning (ML) applications, the demand for explainable artificial intelligence techniques to explain ML models developed for computational chemistry has also emerged. In this study, we present the development of the Boltzmann-weighted cumulative integrated gradients (BCIG) approach for effective explanation of mechanistic insights into ML models trained on high-level quantum mechanical and molecular mechanical (QM/MM) minimum energy pathways. Using the acylation reactions of the Toho-1 β-lactamase and two antibiotics (ampicillin and cefalexin) as the model systems, we show that the BCIG approach could quantitatively attribute the energetic contribution in one system and the relative reactivity of individual steps across different systems to specific chemical processes such as the bond making/breaking and proton transfers. The proposed BCIG contribution attribution method quantifies chemistry-interpretable insights in terms of contributions from each elementary chemical process, which is in agreement with the validating QM/MM calculations and our intuitive mechanistic understandings of the model reactions.

Computational Analysis on the Allostery of Tryptophan Synthase: Relationship between α/β-Ligand Binding and Distal Domain Closure

Tryptophan synthase (TRPS) is a bifunctional enzyme consisting of α and β-subunits and catalyzes the last two steps of l-tryptophan (L-Trp) biosynthesis, namely, cleavage of 3-indole-d-glycerol-3'-phosphate (IGP) into indole and glyceraldehyde-3-phosphate (G3P) in the α-subunit, and a pyridoxal phosphate (PLP)-dependent reaction of indole and l-serine (L-Ser) to produce L-Trp in the β-subunit. Importantly, the IGP binding at the α-subunit affects the β-subunit conformation and its ligand-binding affinity, which, in turn, enhances the enzymatic reaction at the α-subunit. The intersubunit communications in TRPS have been investigated extensively for decades because of the fundamental and pharmaceutical importance, while it is still difficult to answer how TRPS allostery is regulated at the atomic detail. Here, we investigate the allosteric regulation of TRPS by all-atom classical molecular dynamics (MD) simulations and analyze the potential of mean-force (PMF) along conformational changes of the α- and β-subunits. The present simulation has revealed a widely opened conformation of the β-subunit, which provides a pathway for L-Ser to enter into the β-active site. The IGP binding closes the α-subunit and induces a wide opening of the β-subunit, thereby enhancing the binding affinity of L-Ser to the β-subunit. Structural analyses have identified critical hydrogen bonds (HBs) at the interface of the two subunits (αG181-βS178, αP57-βR175, etc.) and HBs between the β-subunit (βT110 - βH115) and a complex of PLP and L-Ser (an α-aminoacrylate intermediate). The former HBs regulate the allosteric, β-subunit opening, whereas the latter HBs are essential for closing the β-subunit in a later step. The proposed mechanism for how the interdomain communication in TRPS is realized with ligand bindings is consistent with the previous experimental data, giving a general idea to interpret the allosteric regulations in multidomain proteins.

Factors That Determine the Variation of Equilibrium and Kinetic Properties of QM/MM Enzyme Simulations: QM Region, Conformation, and Boundary Condition

To analyze the impact of various technical details on the results of quantum mechanical (QM)/molecular mechanical (MM) enzyme simulations, including the QM region size, catechol-O-methyltransferase (COMT) is studied as a model system using an approximate QM/MM method (DFTB3/CHARMM). The results show that key equilibrium and kinetic properties for methyl transfer in COMT exhibit limited variations with respect to the size of the QM region, which ranges from ∼100 to ∼500 atoms in this study. With extensive sampling, local and global structural characteristics of the enzyme are largely conserved across the studied QM regions, while the nature of the transition state (e.g., secondary kinetic isotope effect) and reaction exergonicity are largely maintained. Deviations in the free energy profile with different QM region sizes are similar in magnitude to those observed with changes in other simulation protocols, such as different initial enzyme conformations and boundary conditions. Electronic structural properties, such as the covariance matrix of residual charge fluctuations, appear to exhibit rather long-range correlations, especially when the peptide backbone is included in the QM region; this observation holds when a range-separated DFT approach is used as the QM region, suggesting that delocalization error is unlikely the origin. Overall, the analyses suggest that multiple simulation details determine the results of QM/MM enzyme simulations with comparable contributions.

Retinal Vibrations in Bacteriorhodopsin are Mechanically Harmonic but Electrically Anharmonic: Evidence From Overtone and Combination Bands

Fundamental vibrations of the chromophore in the membrane protein bacteriorhodopsin (BR), a protonated Schiff base retinal, have been studied for decades, both by resonance Raman and by infrared (IR) difference spectroscopy. Such studies started comparing vibrational changes between the initial BR state (all-trans retinal) and the K intermediate (13-cis retinal), being later extended to the rest of intermediates. They contributed to our understanding of the proton-pumping mechanism of BR by exploiting the sensitivity of fundamental vibrational transitions of the retinal to its conformation. Here, we report on new bands in the 2,500 to 1,800 cm⁻¹ region of the K-BR difference FT-IR spectrum. We show that the bands between 2,500 and 2,300 cm⁻¹ originate from overtone and combination transitions from C-C stretches of the retinal. We assigned bands below 2,300 cm⁻¹ to the combination of retinal C-C stretches with methyl rocks and with hydrogen-out-of-plane vibrations. Remarkably, experimental C-C overtone bands appeared at roughly twice the wavenumber of their fundamentals, with anharmonic mechanical constants ≤3.5 cm⁻¹, and in some cases of ∼1 cm⁻¹. Comparison of combination and fundamental bands indicates that most of the mechanical coupling constants are also very small. Despite the mechanical quasi-harmonicity of the C-C stretches, the area of their overtone bands was only ∼50 to ∼100 times smaller than of their fundamental bands. We concluded that electrical anharmonicity, the second mechanism giving intensity to overtone bands, must be particularly high for the retinal C-C stretches. We corroborated the assignments of negative bands in the K-BR difference FT-IR spectrum by ab initio anharmonic vibrational calculations of all-trans retinal in BR using a quantum-mechanics/molecular mechanics approach, reproducing reasonably well the small experimental anharmonic and coupling mechanical constants. Yet, and in spite accounting for both mechanical and electrical anharmonicities, the intensity of overtone C-C transitions was underestimated by a factor of 4–20, indicating room for improvement in state-of-the-art anharmonic vibrational calculations. The relatively intense overtone and combination bands of the retinal might open the possibility to detect retinal conformational changes too subtle to significantly affect fundamental transitions but leaving a footprint in overtone and combination transitions.

Weight average approaches for predicting dynamical properties of biomolecules

Recent advances in atomistic molecular dynamics (MD) simulations of biomolecules allow us to explore their conformational spaces widely, observing large-scale conformational fluctuations or transitions between distinct structures. To reproduce or refine experimental data using MD simulations, structure ensembles, which are characterized by multiple structures and their statistical weights on the rugged free-energy landscapes, are often used. Here, we summarize weight average approaches for various experimental measurements. Weight average approaches are now applied to hybrid quantum mechanics/molecular mechanics MD simulations to predict fast vibrational motions in a protein with a high accuracy for better understanding of molecular functions from atomic structures.

Harder, better, faster, stronger: Large-scale QM and QM/MM for predictive modeling in enzymes and proteins

Computational prediction of enzyme mechanism and protein function requires accurate physics-based models and suitable sampling. We discuss recent advances in large-scale quantum mechanical (QM) modeling of biochemical systems that have reduced the cost of high-accuracy models. Tradeoffs between sampling and accuracy have motivated modeling with molecular mechanics (MM) in a multiscale QM/MM or iterative approach. Limitations to both conventional density-functional theory and classical MM force fields remain for describing noncovalent interactions in comparison to experiment or wavefunction theory. Because predictions of enzyme action (i.e. electrostatics), free energy barriers, and mechanisms are sensitive to the protocol and embedding method in QM/MM, convergence tests and systematic methods for quantifying QM-level interactions are a needed, active area of development.