Quantum-accurate atomistic modeling of enzyme catalysis using a machine learned potential

Introduction

Enzymes are molecular machines that catalyze the chemical reactions underlying life and play a central role in metabolism, signaling, and biotechnology[1, 2, 3]. Scalable mechanistic simulations of enzyme catalysis remain a fundamental challenge because of the requirement for both quantum mechanical (QM) accuracy and the ability to simulate large systems over extended timescales, typically nanoseconds to microseconds [4].

Classical molecular mechanics (MM) force fields rely on a fixed bonding topology and thus cannot model changes in bonding. Reactive force fields permit such changes, but typically require bespoke parameterizations for each system, can miss subtle electronic effects, and tend to reproduce reference reaction energetics only qualitatively [5].

This fundamental trade-off between QM accuracy and biomolecular scale has motivated a long history of approximate modeling [6, 7, 8, 9, 10]. Traditional strategies to reduce complexity include truncated cluster models of the active site [11], which can incur errors from neglecting a large portion of the system, and hybrid QM/MM approaches, which embed a reactive QM region within a larger classical environment [4, 8].

However, such calculations necessitate system-specific knowledge and expert judgment regarding definition of the QM and MM regions, selection of QM functional and MM force fields, and treatment of the QM/MM boundary. This often requires extensive parameter tuning and can risk overfitting to experimental data [8, 12, 13, 14]. The quality of QM/MM predictions is particularly sensitive to choice of the QM region: too small leads to unreliable results, too large becomes computationally prohibitive [14, 15, 16]. Semi-empirical methods [17, 18] can reduce computational cost, typically at the expense of accuracy [19, 20], although new methods are reshaping this trade-off [21]. When solvent plays a role in the reaction, the manually chosen number of QM water molecules adds another parameter that can strongly influence results [22, 23]. QM/MM still remains computationally expensive, with simulations typically limited to tens of picoseconds, [23, 24] making it challenging to observe barrier crossing events.

Machine-learned interatomic potentials (MLIPs) offer an appealing alternative. MLIPs are trained on QM datasets to replicate the performance of QM simulations with orders of magnitude lower computational cost. Large biomolecules represent challenging systems for MLIPs because of their chemically heterogeneous nature, complex many-body effects, and the requirement of stable long-time molecular dynamics (MD) integration. Recent specialized frameworks have shown great promise. Wang et al. introduced AI2BMD, which achieves ab initio protein energies and forces using a fragmentation strategy and treats solvent using a polarizable model [25]. The GEMS and SO3LR frameworks are complementary routes to ab initio accuracy for biomolecular simulations that have been demonstrated on protein folding and solvation dynamics in small systems [26, 27]. ML/MM approaches, in which an MLIP is used for the QM region, have thus been used to alleviate the simulation costs of QM/MM [28, 29, 30], but the challenges associated with ML region specification, MM force field choice, ML/MM boundary treatment, and solvent modeling all remain.

Recent reports indicate that models trained on the Open Molecules 2025 dataset [32] match DFT performance in numerous chemical and biological applications [29, 33, 34, 35]. We show here that eSEN-omol$^{1}$, the base model of UMA [36], can readily scale and be effectively applied to full-chain enzyme catalysis in explicit solvent, treating all atoms at the same level of theory and removing the need for region and boundary selection (Figure 1A). eSEN-omol is a message-passing MLIP trained on OMol25 at the $\omega$B97M-V/def2-TZVPD level of theory [32]. We leverage established enhanced sampling protocols for rare-event simulation to promote transition state crossing and accurately estimate barrier heights in enzymatic reactions [37]. Our integrated distributed inference framework is able to run fully solvated enzyme simulations comprising 32–54k atoms on 8 GPUs. Additional GPUs accelerate simulations and with up to 32 GPUs simulation speeds of up to millions of steps per day on these systems can be realized (Appendix § 6.§ 2) and total simulation times of $\sim1.0\ \mu\mathrm{s}$ are aggregated in this work.

We demonstrate eSEN-omol’s suitability for biological simulations by reproducing sampled cluster DFT single-point energies with high-fidelity, correctly predicting long-range electrostatic behavior in scenarios relevant to proteins in solvent environments, and keeping proteins correctly folded over many nanoseconds of simulation. Figure 1B-E shows a collection of these validations with supporting details provided in Appendix § 6.§ 3. We then study a diverse set of enzymatic systems with a range of different reactive mechanisms to recover free energy surfaces and reaction rate estimates with a single eSEN-omol model, illustrating the accuracy$^{2}$ and transferability of our method. Our results establish a route to QM-accurate simulations of full-chain, explicit solvent enzyme catalysis without system-specific parameterization or tuning (Figure 1A).

Results and Discussion

Claisen Rearrangement in Chorismate Mutase

Chorismate mutase (CM) is a textbook model for enzyme catalysis [38, 39]. It catalyzes the Claisen rearrangement of chorismate to prephenate, a key step in the biosynthesis of tyrosine and phenylalanine. This reaction also proceeds in aqueous solution via the same concerted pericyclic mechanism. Here, we study the Bacillus subtilis CM enzyme (PDB: 3ZQ8 [40]). Experimentally, the enzyme accelerates the reaction by $10^{6}$ to $10^{7}$-fold relative to aqueous solution [41], with the rate enhancement attributed to preferential stabilization of the chair-like transition state through electrostatic interactions involving an arginine residue in the active site [39, 42, 43].

We first examined the intrinsic reactivity in the gas phase using QM calculations at the $\omega$B97M-V/def2-TZVPD level of theory. The QM potential energy profile agrees closely with eSEN-omol predictions, with respective barrier heights of 169.8 kJ mol$^{-1}$ and 169.0 kJ mol$^{-1}$ (Figure S4), providing an initial validation of the pre-trained MLIP. We then introduced explicit water solvent (Figure 2A) and conducted $\sim$46k-atom eSEN-omol simulations at 300 K and 1 bar with On-the-fly Probability Enhanced Sampling (OPES) [44, 45] along a previously reported collective variable (CV) to drive the Claisen rearrangement (Figure 2B, Appendix § 7) [23, 38]. The resulting free-energy surface (FES) in water gives a barrier of $\Delta G^\ddagger_{\mathrm{water}} = (119.81 \pm6.24)\ \mathrm{kJ\,mol^{-1}}$ (Figure 2) in good agreement with experiment (102.55 $\mathrm{kJ\,mol^{-1}}$ [46]) and previous computational studies

(102.09 $\pm$ 2.51 kJ mol$^{-1}$ [22]) (Figure 2D). Finally, we conducted OPES simulations of the reaction within the full-chain, explicitly-solvated enzyme environment and find that the pre-trained eSEN-omol predicts a reaction barrier of $\Delta G^{\ddagger}_{\mathrm{enzyme}}$ = (75.21 $\pm$ 2.72) kJ mol$^{-1}$, in reasonable agreement with the experimentally-reported value of 64.43 kJ mol$^{-1}$ [41] (Figure 2C,D), and corresponding to a $\sim$45 kJ mol$^{-1}$ reduction compared to the barrier in water. The $\sim$10 kJ mol$^{-1}$ discrepancy with respect to experiment is comparable to the errors in the DFT training energies underpinning eSEN-omol. Specifically, we found that the training level of theory $\omega$B97M-V/def2-TZVPD overestimates the CM gas phase barrier height by a similar $\sim$10 kJ mol$^{-1}$ compared to DLPNO-CCSD(T)/def2-QZVPPD. Others have seen similar errors for $\omega$B97M-V barrier heights [47]. Further, the relative barrier reduction $\Delta\Delta G^{\ddagger} = \Delta G^{\ddagger}_{\mathrm{water}} - \Delta G^{\ddagger}_{\mathrm{enzyme}}$, which can mitigate systematic errors, differs by only $\sim$6.5 kJ mol$^{-1}$.

In addition, we computed the rate constants using OPES flooding (OPESf) [49]. In water, we obtained a rate constant of $k_{\mathrm{water}} = (9.01 \pm4.53) \times10^{-9}$ s$^{-1}$ corresponding to an effective activation free energy of $E^{\mathrm{water}}_{a} = (119.69 \pm1.28)$ kJ mol$^{-1}$ (Figure S7). In the enzyme, the calculated rate constant increases to $k_{\mathrm{enzyme}} = (0.26 \pm0.12)$ s$^{-1}$ implying an effective activation free energy of $E^{\mathrm{enzyme}}_{a} = (76.77 \pm1.16)$ kJ mol$^{-1}$. The calculated kinetic activation energies are in excellent agreement with the free energy barriers within $\sim$1.5 kJ mol$^{-1}$, providing strong support for the consistency of the independent thermodynamic and kinetic barrier height estimations. The predicted $k_{\mathrm{enzyme}}/k_{\mathrm{water}} = 3.2 \times10^{7}$ acceleration is in good agreement with the experimentally reported accelerations of $10^{6} - 10^{7}$-fold [41]. These results (Figure 2D) constitute a marked improvement on replicating experimental barriers over previous ML/MM [48] and match the best QM/MM [23] study that uses a density functional fine-tuned to (high level) electronic-structure theory, while improving on the QM/MM result with an untuned functional. While QM/MM methods in general can already achieve high accuracy [22, 50], we envision eSEN-omol’s combination of speed, strong generalization without system specific tuning, and partitioning free setup will accelerate the field of enzymology and more broadly quantum-accurate simulations.

Ester-Bond Hydrolysis in PETase

The PETase enzyme from Ideonella sakaiensis is a promising biocatalyst for depolymerization of poly(ethylene terephthalate) (PET) [51]. PETase employs the canonical catalytic SER-HIS-ASP triad for PET hydrolysis [52, 53] (Figure 3A,B). Structural and biochemical studies suggest that PET hydrolysis proceeds through sequential acylation and deacylation steps [51]. During acylation, the catalytic SER160 attacks the ester bond of PET to form an acyl-enzyme intermediate (AEI), which is subsequently hydrolyzed by a water molecule to regenerate the active enzyme while releasing the product during the deacylation step. A number of QM/MM studies have probed the mechanism and free energy surface associated with this reaction [54, 55, 56, 57, 58, 59], but differences in DFT settings [55, 58], selection of the QM region [58, 60], length of the polymer chains [56], and choice of CVs [54, 58] have led to debate on whether the acylation [24, 55] or deacylation [54, 61] step is rate limiting and whether each reaction proceeds in a single step [54] or through a two-step tetrahedral intermediate [56].

Here, we utilize the solvated PETase structure from Burgin et al. (PDB: 6EQE), as well as the corresponding reaction coordinates (RCs) [54, 62], and PET is represented by a bis(2-hydroxyethyl) terephthalate (BHET) dimer. We conducted $\sim$54k-atom replica exchange umbrella sampling (REUS) [63] to sample the acylation and deacylation pathways (Appendix § 8.§ 3). Transition state region geometries are consistent with the transition state structures reported by Burgin et al. [54].

As illustrated in Figure 3C, we predict an effective free energy barrier of (79.03 $\pm$ 4.2) kJ mol$^{-1}$, derived from the acylation step, to be rate limiting, in good agreement with prior QM/MM studies [24, 54, 55, 56] and experimental measurements for the BHET dimer [51]. We assume, as have others [24, 54, 55, 56], that the acylation and deacylation steps are separable, as diffusion of BHET away from the enzyme effectively renders the acylation step irreversible. Our results are in good agreement with experimental $k_{\mathrm{cat}}$ derived barrier height of (75.19 $\pm$ 0.4) kJ mol$^{-1}$ [51, 64] and that predicted by Burgin et al. of 76.5 kJ mol$^{-1}$ [54]. Compared to Burgin et al., however, our FES exposes some key mechanistic differences. First, our calculations predict acylation, rather than deacylation, to be rate-limiting, as was found by other studies [24, 55]. The highest free energy along the acylation reaction at transition state TS2 is (79.03 $\pm$ 4.2) kJ mol$^{-1}$ (Figure 3D, Table S6), whereas the highest barrier along the deacylation reaction at transition state TS3 is slightly lower at (60.72 $\pm$ 3.0) kJ mol$^{-1}$ (Figure 3E, Table S6). Second, despite using the identical CV and similar umbrella sampling methods, we predict both acylation and deacylation to proceed by a two-step mechanism with well-defined intermediates as opposed to a single concerted step. The acylation intermediate I1 lies $\sim$19 kJ mol$^{-1}$ lower in free energy than TS1 and $\sim$25 kJ mol$^{-1}$ lower than TS2, splitting the nucleophilic addition of the SER160 side chain and the elimination of the BHET leaving group. The deacylation intermediate I2 lies $\sim$27 kJ mol$^{-1}$ lower in free energy than TS3 and $\sim$19 kJ mol$^{-1}$ lower than TS4, splitting the attack of catalytic water and the elimination of SER160 to regenerate the enzyme. Both intermediates show a distinct

Reaction	Feature	Mean RC	$\Delta G$ (kcal mol$^{-1}$)	95% CI (kcal mol$^{-1}$)	$\Delta G$ (kJ mol$^{-1}$)	95% CI (kJ mol$^{-1}$)
Acylation	First transition (TS1)	0.50	17.39	[16.82, 17.97]	72.78	[70.37, 75.19]
Acylation	Intermediate basin 1 (I1)	5.20	12.82	[11.90, 13.74]	53.65	[49.80, 57.49]
Acylation	Second transition (TS2)	8.22	18.89	[17.90, 19.88]	79.03	[74.87, 83.19]
Deacylation	Third transition (TS3)	-1.17	14.51	[13.79, 15.23]	60.72	[57.70, 63.74]
Deacylation	Intermediate basin 2 (I2)	0.85	8.04	[7.67, 8.42]	33.66	[32.07, 35.25]
Deacylation	Fourth transition (TS4)	4.73	12.41	[11.73, 13.08]	51.91	[49.07, 54.74]

tetrahedral geometry centered on the PET ester carbonyl carbon, consistent with the simulations of Guo et al. [56], which employed a higher level of DFT theory (M06-2X/6-31G) than Burgin et al. (DFTB3) [54] who did not find this behavior. This observation is consistent with previous work showing that PBE and DFTB3 tend to underestimate reaction barriers and produce overly smooth free-energy surfaces, whereas methods incorporating exact exchange can better localize charge along bond-breaking and bond-forming coordinates and resolve finer features [56, 58]. Taken together, our results suggest that eSEN-omol, which is trained at the $\omega$B97M-V/def2-TZVPD level of theory, can accurately predict barrier heights and reaction mechanisms at $\sim1000\times$ speedups (Appendix § 6.§ 2) relative to QM/MM methods with a high level of DFT theory.

Enzyme-catalyzed phosphoryl transfer in nucleoside diphosphate kinase

NDPK is a metal-activated phosphotransferase that balances cellular pools of nucleoside triphosphates for signaling, metabolism, and mitochondrial function [65, 66]. It operates through a ping–pong mechanism in which the active-site HIS122 captures the $\gamma$-phosphoryl group of a nucleoside triphosphate (NTP) such as ATP, and donates it to an incoming nucleoside diphosphate (NDP) [67, 68]. We model the NDPK dephosphorylation half-reaction in quantum-accurate simulations with full-chain proteins, explicit solvent, and a metal ion using eSEN-omol. We initialized our eSEN-omol simulations of NDPK from the Dictyostelium discoideum crystal structure (PDB: 1KDN [65]).

We calculated the FES of the phosphoryl transfer reaction using well-tempered metadynamics (WTMetaD) [70] driven in two collective variables: CV1 tracks nucleophilic attack of HIS122 N$\delta$ on the $\gamma$-phosphorus, and CV2 tracks proton transfer from N$\delta$ of HIS122 to the $\gamma$-phosphoryl oxygen (Appendix § 9). The FES suggests that the low-free energy pathway corresponds to a two-step mechanism rather than the concerted proton/phosphoryl transfer previously hypothesized [11] (Figure 4C). In the first step, the proton on the doubly-protonated HIS122 hops to a $\gamma$-phosphoryl oxygen (I→II→III). The protonated $\gamma$-phosphoryl group then transfers through an SN2-like transition state to N$\delta$ of HIS122, completing the dephosphorylation half-reaction (III→IV→V). Once the phosphorylation reaction is completed, the proton previously transferred from HIS122 to $\gamma$-phosphoryl is free to hop back and forth between the $\beta$-phosphoryl on ADP and the phosphoryl group on HIS122 (V→VI) with a barrier of just $\sim7.5$ kJ mol$^{-1}$. The transition state (VII) of the concerted proton/phosphoryl transfer mechanism (I→VII→V) stands $\sim100$ kJ mol$^{-1}$ above the reactant, which is substantially higher than the corresponding value of $\sim70$ kJ mol$^{-1}$ for the transition state (IV) of the rate-limiting second step of the two-step pathway (III→IV→V). We note that, as is the case for all low-dimensional FES projections, barrier heights depend on the choice of CVs and care should be taken not to overinterpret reaction paths and activation energies. Nevertheless, the CVs here were chosen from physical understanding of the mechanism, resolve both the two-step and concerted pathways within the same FES, and identify the two-step pathway as lower in free energy. As we show below, we observe good agreement of the predicted barrier height for the putative rate limiting step with an OPE$\mathrm{S}_f$ rate calculation.

As an independent verification of the free energy barriers calculated from 2D WTMetaD, we conducted 1D umbrella sampling (US) [71] of the two elementary steps in the proposed two-step pathway. US calculations in CV2 between states I and III expose proton hopping as a largely uphill transition through the transition state II with an elementary free energy barrier of $\Delta G_1^{\ddagger} = (27.49 \pm1.05)$ kJ mol$^{-1}$ (Figure 4D). An analogous US calculation in CV1 between states III and V predicts phosphoryl transfer to be a single-barrier two-state mechanism with an elementary free energy barrier of $\Delta G_2^{\ddagger} = (73.05 \pm1.86)$ kJ mol$^{-1}$ (Figure 4E). The barriers computed by 1D US are in agreement with those computed by marginalization over the reaction channels in the 2D WTMetaD FES (Figure S18). The overall free energy change for the dephosphorylation half-reaction between states I and V is approximately $\Delta G = 15$ kJ mol$^{-1}$.

We applied OPE$\mathrm{S}_f$ to the rate limiting phosphoryl transfer step (Figure 4E, Section A.4) [49, 72] and obtained a predicted rate constant of $k_{\mathrm{OPE S_f}} = (0.043 \pm0.019)$ s$^{-1}$ with a corresponding apparent effective activation free energy of $E_a^{\mathrm{OPE S_f}} = (81.32 \pm1.13)$ kJ mol$^{-1}$ (Figure S19). The effective activation energy agrees with the free energy barrier computed by US to within $\sim10$ kJ mol$^{-1}$. In Figure 4F we compare our US (73 kJ mol$^{-1}$) and OPE$\mathrm{S}_f$ (80.3 kJ mol$^{-1}$) barriers to the value implied by the experimental $k_{\mathrm{cat}}$ (55.6 kJ mol$^{-1}$ [69]). Because $k_{\mathrm{cat}}$ reports on the full cycle, which is diffusion-limited, this value provides an upper bound on the chemical barrier. The US estimate shows the best agreement with experiment, exceeding it by $\sim17$ kJ mol$^{-1}$, representing a roughly five-fold error reduction from the previous semi-empirical AM1

cluster-model (148.4 kJ mol$^{-1}$ [11]).

Conclusion

The ability to perform quantum-mechanically accurate simulations of full-chain, explicit-solvent enzyme catalysis without system specific setup or tuning has the potential to transform mechanistic understanding and rational design of enzyme catalysis. In this work, we demonstrate that pre-trained MLIPs provide a viable pathway to this goal. eSEN-omol, trained on OMol25, simultaneously offers QM-accurate calculations at the $\omega$B97M-V/def2-TZVPD level of theory and $\sim$1000$\times$ speedups over equivalent QM/MM methods (Appendix A.2), as well as treatment of full-chain enzymes in explicit solvent at a unified level of theory without any system-specific setup and tuning. We couple the eSEN-omol simulations to off-the-shelf rare event enhanced sampling techniques – REUS [63], WTMetaD [70], OPES [45], and OPESf [49] – to efficiently estimate free energy surfaces, barriers, and rates. In an application to chorismate mutase, we accurately predicted effective barriers and rate accelerations consistent with experimental measurements. In PETase, our rate-limiting barriers are in good agreement with experiment and prior QM/MM calculations and we resolve previously reported metastable intermediates. In the NDPK dephosphorylation half-reaction, we predicted a two-step reaction mechanism mediated by interstitial solvent molecules and a Mg$^{2+}$ cation in which proton transfer is followed by a rate-limiting phosphoryl transfer, yielding a reaction barrier with an approximately five-fold lower deviation from experiment than a previous semi-empirical QM study.

While transferable MLIPs have demonstrated great promise across broad areas of chemistry [32] and we demonstrated here their promise in biomolecular simulations of enzymatic catalysis, these methods are still relatively new and have certain current limitations and open questions. MLIPs inherit any shortcomings that exist in the training data, including systematic errors associated with the level of DFT theory and statistical errors associated with the diversity of the training set. MLIPs are also generally untested on systems which may require careful treatment of significant charge transfer or multiple electronic or spin states. Unlike classical force fields, long-range electrostatics and induced polarization are implicitly encoded in a graph neural network and not directly accessible for modulation. Further, while MLIPs are much faster than QM calculations, they remain significantly slower than classical potentials, so enhanced sampling techniques with their attendant complexities are still required for MLIP-driven simulations of enzyme catalysis.

The applications in this work have demonstrated that MLIPs can open new avenues for scientific discovery in mechanistic understanding and computational design of enzymes going beyond the prediction of stable structures [73, 74, 75, 76] and conformational ensembles [77, 78, 79, 80] to fully embrace the dynamical complexity of enzyme catalysis. In the NDPK application, for example, eSEN-omol enabled the prediction of complex reaction mechanisms that would otherwise remain difficult to model or suffer from strong biases associated with system-specific setup. As these methods continue to mature, they promise to transform our ability not only to understand complex biomolecular processes, but also to systematically design systems with targeted functions.

Acknowledgements

We thank Juno Nam for fruitful conversations and critical feedback on an early draft of the manuscript. We additionally thank Bingqing Cheng for providing dipole correlation data. Calculations were completed in part with resources provided by the University of Chicago Research Computing Center. We gratefully acknowledge computing time on the University of Chicago high-performance GPU-based cyberinfrastructure supported by the National Science Foundation under Grant No. DMR-1828629.

Funding: SKA acknowledges support from the Eric and Wendy Schmidt AI in Science Fellowship, a Schmidt Sciences program.

Author contributions:

Conceptualization: ASZ, CLZ, DSL, ALF;

Data curation: ASZ, MG, AS, SD, SKA, MD, BMW, DSL;

Formal analysis: ASZ, MG, AS, SD, SKA; Investigation: ASZ, MG, AS, SD, SKA, MD, BKM, CLZ, BMW, ZWU, DSL, ALF;

Methodology: ASZ, MG, AS, SD, SKA, ALF, BMW;

Project administration: MG, ZWU, DSL, ALF;

Resources: CLZ, ZWU;

Software: MG, MD, BMW;

Supervision: CLZ, ZWU, ALF;

Validation: ASZ, MG, AS, SD, SKA, BMW, DSL;

Visualization: ASZ, MG, AS, SD, SKA, BMW;

Writing – original draft: ASZ, MG, AS, SD, SKA, MD, BMW, DSL, ALF;

Writing – review & editing: MG, BMW, DSL, ALF.

Competing interests: ALF is a co-founder and consultant of Evozyne, Inc. and a co-author of US Patent Applications 16/887,710 and 17/642,582, US Provisional Patent Applications 62/853,919, 62/900,420, 63/314,898, 63/479,378, 63/479,378, 63/521,617, 63/510,130, 63/669,836, and 63/987,554 and International Patent Applications PCT/US2020/035206, PCT/US2020/050466, PCT/US24/10805, PCT/US24/34369, and PC-T/US25/35833.

Data, code, and materials availability: All data are available in the manuscript or the supplementary materials. The eSEN models are available at https://huggingface.co/facebook/OMol25. Code for running models is available at https://github.com/facebookresearch/fairchem.

Supplementary Materials

Materials and Methods

Supplementary Text, Figures, Tables, and References

[1] T. Hunter, Protein kinases and phosphatases: The Yin and Yang of protein phosphorylation and signaling. Cell 80 (2), 225–236 (1995), doi:10.1016/0092-8674(95)90405-0.

[2] U. T. Bornscheuer, et al., Engineering the third wave of biocatalysis. Nature 485 (7397), 185–194 (2012), doi:10.1038/nature11117.

[3] P. K. Robinson, Enzymes: principles and biotechnological applications. Essays in Biochemistry 59, 1–41 (2015), doi:10.1042/bse0590001.

[4] A. Warshel, Multiscale Modeling of Biological Functions: From Enzymes to Molecular Machines (Nobel Lecture). Angewandte Chemie International Edition 53 (38), 10020–10031 (2014), doi:https://doi.org/10.1002/anie.201403689.

[5] T. Trnka, I. Tvaroška, J. Koča, Automated Training of ReaxFF Reactive Force Fields for Energetics of Enzymatic Reactions. Journal of Chemical Theory and Computation 14 (1), 291–302 (2018), doi:10.1021/acs.jctc.7b00870.

[6] X. Lu, J. Chen, J. Huang, The continuous evolution of biomolecular force fields. Structure 33 (7), 1138–1149 (2025), doi:10.1016/j.str.2025.05.013.

[7] S. J. Marrink, et al., Two decades of Martini: Better beads, broader scope. WIREs Computational Molecular Science 13 (1), e1620 (2023), doi:10.1002/wcms.1620.

[8] H. M. Senn, W. Thiel, QM/MM Methods for Biomolecular Systems. Angewandte Chemie International Edition 48 (7), 1198–1229 (2009), doi:10.1002/anie.200802019.

Appendix

Materials and Methods

Molecular Dynamics Simulations Using eSEN-omol

Molecular modeling calculations were conducted using the ASE suite [81] with the eSEN-omol MLIP through Meta’s fairchem package [82]. Enhanced sampling simulations were performed with PLUMED 2.11.0 enhanced sampling libraries [83].

Distributed Inference

In order to perform enhanced sampling/molecular dynamics of these large enzymatic systems, we need to be able to run simulations on nanosecond/day timescales with $O(100\mathrm{k})$ atoms with a highly accurate model for molecules. There are many MLIPs that can run simulations on this number of atoms but they are either not accurate enough for this application, do not support parallelism [84], rely on strictly local interactions [85, 86], or are limited to inter-node only parallelism [87]. Towards this end, we designed our domain decomposed parallelism algorithms to scale eSEN-omol to millions of steps/day (equivalent to ns/day for 1 fs timesteps) speeds for 100k atom systems with any simulator engine including ASE and LAMMPS. The code is freely available on GitHub at https://github.com/facebookresearch/fairchem/. Execution speeds with ASE are reported in Table S1.

System	Atoms	1 GPU	8 GPU	32 GPU
Water (1000 mol.)	3,000	0.86	3.17	3.27
NDPK	32,339	oom	0.60	1.53
Chorismate mutase	46,323	oom	0.38	1.06
PETase (acylation)	53,816	oom	0.35	0.97

In this paper, we reference 1000x speedup relative to high level QM/MM methods (for example: M06-2x, $\omega$B97X-3c). This is an approximate argument based on what is frequently reported in QM/MM studies where methods with $\omega$B97X-3c report 3000-6000 steps per day[23] and M06-2x report total simulation times of approximately tens of ps[56] for example.

Note. Single-GPU runs for the solvated enzymes exhaust device memory, as expected at these system sizes.

Suitability of eSEN-omol for Large Scale Biological Simulations

DFT Accuracy

We demonstrate that eSEN-omol matches DFT to within a few kJ/mol accuracy in a 25 Å sampled QM region for configurations harvested from an enzyme catalyzed reactive transition from chorismate to prephenate by the chorismate mutase enzyme under a 1.42 ns, 300 K OPES metadynamics trajectory (Figure 1B, cf. § 2.1). Defining the collective variable $\Delta d = d(\mathrm{C16}-\mathrm{C3}) - d(\mathrm{O11}-\mathrm{C1})$, where the two terms describe the forming of C-C and breaking of C-O bonds in the CHB substrate, fifty target values were uniformly spaced between $\Delta d = (-0.060)$-0.156 nm. For each target, the saved trajectory frame having the nearest CV value was selected and each frame was permitted to be used only once. Frames were available every 500 MD steps (0.25 ps). The selected configurations spanned 17.0-1400.75 ps, and the maximum deviation from a target CV value was (5.33 $\times10^{-4}$) nm. A fixed-composition active-site cluster was constructed from each configuration. Every cluster contained the complete CHB substrate; protein residues (denoted by chain name and residue) B57–B60, B63, B73–B75, C7–C9, C78, C90, and C108; and one conserved water molecule (SPC 8334). The catalytic residues were selected following previous QM/MM descriptions of the inter-subunit active site [88, 89], with additional persistent contact and packing residues included subject to a 350-atom computational limit. Pro B58 and Gly C8 were retained to avoid introducing adjacent, overlapping peptide caps. The water molecule was the most persistent water near CHB; its oxygen-to-substrate-heavy-atom distance ranged from 2.55 to 5.00 Å among the selected configurations. The retained protein regions formed seven peptide fragments. Each fragment was terminated with an acetyl group at its N terminus and an N-methylamide group at its C terminus. Cap heavy-atom positions were inherited from the instantaneous geometry of the neighboring backbone, while cap hydrogens were placed using idealized bond lengths and tetrahedral geometry. All other coordinates were transferred directly from the trajectory without geometry optimization. Each resulting cluster contained 343 atoms, had a formal charge of +1, and was treated as a singlet. The atom composition and ordering were identical for all configurations. As validation, the CV was recomputed directly from each trajectory frame and agreed with the corresponding COLVAR value within $(2 \times10^{-6})$ nm. Hydrogen connectivity was checked at every selected geometry, and no interatomic separation below 0.70 Å was permitted. These CV-stratified snapshots were randomly subsampled to choose 24 structures. Note that orthogonal environmental coordinates are uncontrolled, so absolute cluster energies should not be interpreted directly as a one-dimensional reaction profile.

Simulation Stability

We also demonstrate that all-atom simulations of full-chain chorismate mutase, PETase, and NDPK proteins in explicit solvent and containing a significant number of charged residues, charged substrates, and ions using eSEN-omol are stable and remain properly folded ($\mathrm{C}_{\alpha}$ RMSD $< 2$ Å) over multi-nanosecond time scales (Figure 1C, S2, S8, and S14). While the time scales of these unbiased calculations are insufficiently long to observe meaningful conformational changes in large biomolecular systems, the stability of MLIPs over nanosecond time scales is not necessarily given or even assessed in many applications [90].

Evidence for Correct Long Range Behavior

Highly charged peptide. A 2 ns simulation of the highly charged 1FME peptide [91] comprising 14 charged residues (out of 28 total) and numerous long-range salt bridges produces a close match to experimental NMR structures (Figure 1C) and classical force field results (AMBER FF19SB). We solvated 1FME in water (12k atoms in total), minimized the energy, equilibrated 100 ps of NVT at 300K, 100 ps of NPT at 1 atm and then ran unbiased simulations for 2 ns. The $\mathrm{C}_{\alpha}$ RMSD is averaged over the run.

Dipole density correlations. We show that our model produces the correct long range dipole correlations in water (Figure 1D). Replicating the calculations from Ref. [31], we perform NVT simulations of a box with 512 atoms over 2 ns and show the long range dipole correlation function $\langle\tilde{\mathbf{m}}^{*}(\mathbf{k})\tilde{\mathbf{m}}(\mathbf{k})\rangle$, where $\tilde{\mathbf{m}}$ is the Fourier transform of the molecular dipole density, agrees with the long-range models with explicit Ewald summation out to the $2\pi/L$ k-space limit. Specifically, we compute the molecular dipole density of water,

where $\boldsymbol{\mu}_j$ is the dipole of water molecule $j$ and $\mathbf{r}_j$ is the position of the water oxygen atom. For each allowed wavevector in reciprocal space, $\mathbf{k} = \frac{2\pi}{L}(n_x,n_y,n_z)$, we take the dipole moment projection,

then we take the Fourier transform to compute the dipole density in k-space,

where the reported value is averaged over trajectory frames. k-vectors of the same magnitude are averaged together. In order to draw quantitative comparisons with [31], we use the same box size with 512 atoms and a length of 24.86 Å which limits the smallest k computed to be 0.25 Å$^{-1}$. We then overlaid Cheng’s results on ours, shifting the spectra by a constant offset along the y-axis to align the first peak.

PMF of ions in solution. We show that eSEN-omol correctly captures the long-range behavior of ions in solution. The high dielectric constant of water ($\epsilon\approx78$ at room temperature) means that charge interactions are screened rapidly over short distances and the potential of mean force (PMF) of a NaCl ion pair closely agrees with force fields with explicit long range electrostatics terms (Figure 1E). In a low dielectric solvent like hexane, we show that eSEN-omol captures electrostatic interactions in the NaCl ion pair out to $>20$ Å ionic separation (at which point the potential is effectively 0 kJ/mol) in excellent agreement with both the continuum limit and classical simulations using the CHARMM-Drude force field (Figure 1F).

The classical reference system contained one Na$^{+}$, one Cl$^{-}$, and 900 water molecules in a fixed 30 Å cubic cell at 300 K. Three classical models were considered: Joung–Cheatham ions with rigid TIP3P water, polarizable AMOEBA2018 water and ions, and CHARMM-Drude ions with SWM4-NDP water. Umbrella sampling employed 16 windows spanning Na–Cl separations of 2.0–9.0 Å, with a force constant of 20 kJ mol$^{-1}$ Å$^{-2}$. Following 100 ps of restrained equilibration, production lengths were 2 ns, 1 ns, and 0.5 ns per window for JC/TIP3P, AMOEBA2018, and CHARMM-Drude, respectively.

In low-dielectric $n$-hexane, the reference system contained one NaCl pair and 691 hexane molecules in a fixed 53.25 Å cubic cell at 300 K. The classical calculation used the polarizable CHARMM-Drude 2023 force field with matching ion and all-atom hexane parameters, PME electrostatics, a 12 Å nonbonded cutoff, and a 10 Å switching distance. Umbrella windows covered separations from 2.5 to 20 Å, using progressively weaker restraints and wider spacing at larger separations; each window included 25 ps of restrained equilibration and 300 ps of production. The eSEN-omol and CHARMM-Drude PMFs show similar distance dependence across the sampled range (Figure 2E). Moreover, fitting eSEN-omol’s PMF, to an $A/r + C$ functional form over 8–18 Å gave, gave an effective dielectric constant of approximately 1.87, consistent with the experimental continuum value of 1.88 for hexane.

All PMFs were reconstructed using BayesWHAM [92], corrected for the radial Jacobian by adding $2k_{\mathrm{B}}T\ln r$, and referenced to zero over 8–9 Å in water and 19–20 Å in hexane. Uncertainties were estimated using time-block analyses or moving-block bootstrap reconstructions. The zero value near 20 Å is therefore a chosen free-energy reference, rather than evidence that the physical ion–ion interaction has vanished completely. The short-range CHARMM-Drude water PMF below 3.25 Å was treated as qualitative because of activity of the model’s Drude-particle hard wall.

In Appendices § 10 and § 11, we provide additional assessments of the performance of eSEN-omol in predicting the structure and dynamics of liquid water and the folding of a small peptide.

Taken together, these results provide support that eSEN-omol is well-suited to protein-scale, condensed phase calculations in bulk water. Despite not having an explicit Coulombic energy term commonly found in classical force fields, it accurately captures long range electrostatic interactions over its effective receptive field of $r_{\mathrm{eff}} = 24$ Å that is sufficient to capture highly screened charge interactions in water solvent.

Rate Estimation with Eyring’s Transition State Theory and OPES Flooding

To quantify reaction kinetics, we employed transition-state theory (TST) and OPES flooding (OPESf). Within Eyring’s TST, the rate constant associated with an activation free-energy barrier $\Delta G^{\ddagger}$ is given by,

where $\kappa$ is the transmission coefficient, here assumed to be $\kappa= 1$, $k_{\mathrm{B}}$ is Boltzmann’s constant, $h$ is Planck’s constant, and $T$ is the temperature. The equilibrium free-energy barriers $\Delta G^{\ddagger}$ obtained from enhanced sampling calculations neglect dynamical recrossings of the transition-state dividing surface. Accordingly, we computed reaction rates using OPESf, which directly samples barrier-crossing events and thereby incorporates dynamical effects beyond the equilibrium TST estimate.

OPESf accelerates rare transitions by progressively flooding the reactant basin while leaving the transition-state region unbiased [93]. A history-dependent bias $V_b(s,t)$ is deposited along a collective variable $s$, while an excluded region prevents bias deposition beyond a system-specific threshold $s_{\mathrm{exc}}$. The maximum bias is controlled by the energy cutoff $\Delta E$. The parameters $\Delta E$ and $s_{\mathrm{exc}}$ were selected so as to reduce the effective barrier within the reactant basin while preserving unbiased dynamics in the transition-state region.

Provided that the transition-state region remains unaffected by the bias, the biased first-passage time $t_f$ can be rescaled to the corresponding unbiased transition time $t^*$ according to,

where $\beta= 1/k_{\mathrm{B}}T$, $U$ is the unbiased potential energy, and the ensemble average defines the acceleration factor associated with the flooding bias [93, 94, 95, 96, 97]. For a rare-event process, the distribution of unbiased transition times is expected to follow a Poisson process with cumulative distribution function,

where $\tau$ is the characteristic transition time. The reaction rate was obtained as,

by fitting the distribution of $t^*$ values obtained from an ensemble of independent OPESf trajectories. The consistency of the transition-time distribution with a Poisson process was assessed using a two-sample Kolmogorov–Smirnov test [97, 98, 99].

Uncertainties in the OPESf transition times and rate constants were estimated following the procedure described in Ref. [72]. The 95% confidence interval $(\tau_l,\tau_u)$ for the characteristic transition time was obtained using Kaminsky’s method,

where $t_i$ is the rescaled transition time for the $i$th trajectory and $\chi^2_{2n}(\alpha)$ denotes the $\alpha$-quantile of the chi-squared distribution with $2n$ degrees of freedom. A symmetric uncertainty was reported as,

Since $k = 1/\tau$, the corresponding confidence limits on the rate are,

with the reported symmetric uncertainty,

The OPESf rate constants were converted to effective activation barriers by inverting the Eyring equation,

The confidence limits on the effective barrier were obtained as,

A symmetric uncertainty in the effective activation barrier was then reported as,

such that the OPESf-derived effective barrier is expressed as $\Delta G_{\mathrm{eff}}^{\ddagger} \pm\Delta\left(\Delta G_{\mathrm{eff}}^{\ddagger}\right)$. The OPESf parameters employed for the different systems are summarized in Table S2.

System	Collective variable	$\Delta E$ (kJ mol$^{-1}$)	Product region	EXCLUDED_REGION
CM (Aqueous)	CV	120	CV $> 0.15$ nm	CV $< -0.15$ nm
CM (Enzyme)	CV	60	CV $> 0.15$ nm	CV $< -0.15$ nm
NDPK	CV$_1$	60	CV$_1 > 0.4$ nm	CV$_1 < -0.32$ nm

Chorismate Mutase

Simulation Setup and Details

Initial structures for the full enzyme were obtained from Ray et al. [100], with the protein conformation based on PDB 3ZO8 [40]. For the corresponding aqueous simulations, the system consisted only of the reactive chorismate (CHB) species, with the protein removed and the solute fully solvated in water. In all chorismate mutase simulations, we employed ASE [81] together with eSEN-omol [36]. Prior to the OPES simulations, each system was first energy minimized and equilibrated. For energy minimization, we used a force tolerance of $0.05$ eV/Å. We performed a two-stage equilibration, first in the NVT ensemble at $T = 300$ K, followed by equilibration in the NPT ensemble at $T = 300$ K and $P = 1$ bar. We employed a Langevin thermostat during the equilibration and final scale simulations and an isotropic MTKNPT barostat during the NPT equilibration.

Enhanced Sampling Calculations

Following the reaction coordinate used previously for OPES flooding calculations of chorismate mutase [23], the Claisen rearrangement of chorismate to prephenate was described using a single CV defined as $\mathrm{CV} = d_1 - d_2$, where $d_1$ is the distance associated with formation of the new C–C bond and $d_2$ is the distance associated with cleavage of the C–O bond (Figure 2). Progression from reactant to product is captured by the concerted decrease in $d_1$ and increase in $d_2$. OPES simulations were performed at $300$ K using a bias deposition pace of 500 steps and a bias barrier of 75 kJ mol$^{-1}$, with the adaptive kernel width updated every 1000 steps. Harmonic walls were applied along the reaction coordinate to restrict sampling to the region relevant to the forward reaction, with a lower wall at $\mathrm{CV} = (-0.05)$ nm and an upper wall at $\mathrm{CV} = 0.27$ nm. Both walls used a force constant of $5.0 \times10^{4}$ kJ mol$^{-1}$ nm$^{-2}$. The same OPES setup, but without the lower wall, was used for the corresponding aqueous-phase reaction, in which the reactive substrate was simulated in explicit solvent without the enzyme.

Convergence of the OPES simulation is assessed in Figure S3. The distances $d_1$, $d_2$, and the CV undergo multiple spontaneous transitions between the reactant and product basins throughout the 1.18 ns trajectory (Figure S3–Figure S3), satisfying the multiple-crossing criterion for sampling convergence [45, 101]. The number of deposited OPES kernels plateaus at $\sim$111 (Figure S3), indicating that the adaptive bias has ceased to grow and the underlying free-energy landscape is well sampled. The lower and upper harmonic wall biases (LB and UB; Figure S3–Figure S3) are activated only transiently, confirming that the OPES bias alone drives barrier crossing. and the walls serve exclusively as boundary safeguards. Free-energy surfaces reweighted from five cumulative fractions of the trajectory (20%–100%) converge to a consistent barrier of $\sim109$ kJ mol$^{-1}$, with all estimates within $\sim3$ kJ mol$^{-1}$ of the full-data reference (Figure S3).

Going beyond equilibrium free energies, we next computed rate constants using OPES flooding (OPESf) [93]. The reweighted first-passage times of ensembles of 15 and 18 independent OPESf trajectories collected, respectively, for the aqueous and enzyme-catalyzed reactions follow Poisson statistics, from which we extracted estimates of the unbiased reaction rates (Section A.4, Figure S7) [49, 72]. In water, we computed a rate constant of $k_{\mathrm{water}} = (9.01 \pm4.53) \times10^{-9}\ \mathrm{s}^{-1}$ corresponding to an effective activation free energy of $E_a^{\mathrm{water}} = (119.69 \pm1.28)$ kJ mol$^{-1}$ (Figure S7). In the enzyme, the calculated rate constant increases to $k_{\mathrm{enzyme}} = (0.26 \pm0.12)\ \mathrm{s}^{-1}$ with a corresponding effective activation free energy of $E_a^{\mathrm{enzyme}} = (76.77 \pm1.16)$ kJ mol$^{-1}$. The calculated kinetic activation energies are in excellent agreement with the free energy barriers within $\sim1.5$ kJ mol$^{-1}$, providing strong support for the consistency of the independent thermodynamic and kinetic barrier height estimations. Our kinetic calculations predict a $k_{\mathrm{enzyme}}/k_{\mathrm{water}} = 3.2 \times10^7$ acceleration of the enzyme-catalyzed reaction relative to aqueous solution, in good agreement with the experimentally reported accelerations of $10^6$ to $10^7$-fold [41].

Structural Model Variation and Comparison to Li and Chan Results

We note that when we attempted to dock chorismate into the PDB: 2CHT [102] crystal structure, we initially obtained an improperly docked system in which the reaction center is largely solvent exposed and the substrate is only stabilized by a single hydrogen bond (Figure S6). Our eSEN-omol OPESf calculations for this structure yield a similarly elevated free energy barrier of 108 kJ mol$^{-1}$ (Figure S5) compared to Li and Chan’s results on their 1HB structure [23]. In contrast, the higher-resolution PDB: 3ZO8 [40] structure resulted in a more stably docked substrate, featuring a two-hydrogen-bonded substrate conformation and a substantially lower activation barrier, directly supporting the original observations from Li and Chan [23]. This result also highlights that while studies of reaction dynamics are now readily feasible with MLIPs, an accurate initial structural model is still a vital first step.

We also note that in Figure 2, we only compared to Li and Chan’s values on $\omega$B97X-3c to keep the level of theory and functional comparable to eSEN-omol. Li and Chan showed that refining $\omega$B97X-3c parameters on LNO-CCSD(T) on chorismate single points indeed improved their results. In Section B.4, we find that computing barriers with DLPNO-CCSD(T) also yields similar improvements, suggesting that future models should be trained on CCSD(T) data when they become available.

Ab initio validation of gas phase barrier height prediction

We have quantified the error in barrier height predicted by eSEN-omol, relative to its own reference level of theory ($\omega$B97M-V/def2-TZVPD). The results are shown in Figure S4 and display excellent agreement. To quantify the error in the barrier height relative to reference quantum chemistry (gold-standard) estimates, we have also recomputed the barrier with DLPNO-CCSD(T)/def2-QZVPPD using the geometry from $\omega$B97M-V/def2-TZVPD; the results are shown in Table S3. The ab initio calculations were performed with ORCA version 6.0 and used $\omega$B97M-V/def2-TZVPD settings identical to the OMol25 dataset [32]. Concretely, these calculations use the $\omega$B97M-V functional in conjunction with the def2-TZVPD basis set, the resolution of identity (RI) approximation and tight SCF convergence settings. In addition to the SCF settings employed in the $\omega$B97M-V/def2-TZVPD calculations we use the TightPNO setting for DLPNO. The reactant and product structures were optimized at the $\omega$B97M-V/def2-TZVPD level; subsequently the TS geometry was obtained using the growing string method and optimized with ORCA’s TS Opt. The transition state structure was verified at $\omega$B97M-V/def2-TZVPD level to contain one imaginary frequency.

Level of theory	Barrier height (kj/mol)	Reaction energy (kj/mol)
$\omega$B97M-V/def2-TZVPD	169.8	-90.5
eSEN-omol single-points	170.5	-89.0
eSEN-omol relaxations	169.0	-89.0
DLPNO-CCSD(T)/def2-TZVPPD	156.2	-91.0
DLPNO-CCSD(T)/def2-QZVPPD	159.4	-89.8

PETase

Simulation Setup and Details

Initial structures for both the acylation and deacylation simulations were obtained from [54]. The deacylation simulations were initiated from the same acyl-enzyme intermediate (AEI) structure used in [54]. In all the PETase simulations, we employed ASE [81] together with eSEN-omol [36]. Prior to the WTMetaD simulations, each system was first energy minimized and equilibrated. For energy minimization, we used a force tolerance of $0.05\ \mathrm{eV}/\mathring{\mathrm{A}}$. We performed a two-stage equilibration to gradually bring the system to the target temperature and density. In the first stage, the system was equilibrated in the NVT ensemble at $T = 310\ \mathrm{K}$. Subsequently, the system was equilibrated in the NPT ensemble at $T = 310\ \mathrm{K}$ and $P = 1\ \mathrm{bar}$. We employed the Langevin thermostat and IsotropicMTKNPT barostat during the equilibration and final scale simulations.

Reaction Coordinates

For all simulations in the PETase section, we used identical Reaction Coordinates (RCs) reported by [54] for both the acylation and deacylation reactions. The acylation reaction coordinate ($\mathrm{RC}_{\mathrm{acylation}}$) is a linear combination of three collective variables ($\mathrm{CV1a}$, $\mathrm{CV2a}$, $\mathrm{CV3a}$). $\mathrm{CV1a}$ and $\mathrm{CV2a}$ correspond to the distance differences $d_{\mathrm{SER160:H-PET:O}} - d_{\mathrm{SER160:O-PET:C}}$ and $d_{\mathrm{SER160:O-PET:C}} - d_{\mathrm{PET:C-PET:O}}$, respectively. $\mathrm{CV3a}$ corresponds to the angle formed by SER160 H–PET ether O–BHET C. Similarly, the deacylation reaction coordinate ($\mathrm{RC}_{\mathrm{deacylation}}$) is a linear combination of three collective variables ($\mathrm{CV1d}$, $\mathrm{CV2d}$, $\mathrm{CV3d}$). $\mathrm{CV1d}$, $\mathrm{CV2d}$, and $\mathrm{CV3d}$ correspond to the distance differences $d_{\mathrm{SER160:O-Wat:H}} - d_{\mathrm{Wat:H-Wat:O}}$, $d_{\mathrm{SER160:O-MHET:C}} - d_{\mathrm{MHET:C-Wat:O}}$, and $d_{\mathrm{MHET:C-Wat:O}} - d_{\mathrm{Wat:O-Wat:H}}$, respectively. These reaction coordinates collectively describe the proton-transfer, nucleophilic attack, and bond-breaking/bond-forming events during the acylation and deacylation reactions. The RCs are:

Replica Exchange Umbrella Sampling

For mapping the free energy surface of the full reaction during both acylation and deacylation, we used replica exchange umbrella sampling (REUS) [63]. In this variant of umbrella sampling (US) [71], exchange of configurations between neighboring windows is allowed through a Metropolis criterion to improve convergence and sampling of orthogonal degrees of freedom. In particular, for a given window $i$, assuming a 1D CV $q(x)$, the total energy of the system in REUS is described with $H_i(x,p) = U(x) + K(p) + w_i(q(x))$ where $w_i(q) = \frac{1}{2}k_i(q-c_i)^2$ represents a harmonic bias potential added to restrain and sample the system near the window center $c_i$ with a force constant of $k_i$. $H_i(x,p)$, $U(x)$, and $K(p)$ represent the Hamiltonian, potential energy, and kinetic energy of the system, respectively. $q$ and $p$ refer to the configuration and momentum of the system, respectively. The configurations and momenta from neighboring windows $i$ and $j$ in REUS are exchanged every 1000 simulation steps by using the Metropolis criterion with acceptance probability $P_{\mathrm{acc}} = \min[1,\exp(\beta\delta)]$ where $\delta= w_i(q_i) + w_j(q_j) - w_i(q_j) - w_j(q_i)$ and $\beta= \frac{1}{k_{\mathrm{B}}T}$. When $\delta\ge0$, the swap lowers the overall energy and is always accepted.

We obtain the starting configurations for each window by computing an approximate profile of the full free energy surface using OPES MetaD [44] simulation by patching ASE with PLUMED [83]. In Figure S9 and Figure S10, we show the resulting histograms of $RC_{\mathrm{acylation}}$ and $RC_{\mathrm{deacylation}}$ in each window during the acylation and deacylation REUS runs, respectively. We adjusted window centers ($c_i$) and the value of the spring force constant ($k_i$) as needed to restrain the system in a given window and improve the overlap of the biased histograms of the neighboring windows in both $RC_{\mathrm{acylation}}$ and $RC_{\mathrm{deacylation}}$ as shown in Table S4 and Table S5, respectively. These tables also show the length of the simulation performed in each window. We compute the free energy profiles from the resulting biased distributions of $RC_{\mathrm{acylation}}$ and $RC_{\mathrm{deacylation}}$ in each window using BayesWHAM [92, 104]. We detail the key features on the free energy surfaces in Table S6. The total stable sampling time was $\sim400$ ps per window, totaling $\sim50$ ns of simulation time. The samples were divided into 6 even blocks to obtain measurement errors.

ID	Center	$k$	ID	Center	$k$	ID	Center	$k$
0	-10	20	71	-0.15	75	35	4.5	20
1	-9.5	20	72	-0.1	150	36	5	20
2	-9	20	73	-0.05	250	37	5.5	20
3	-8.5	20	69	0	300	38	6	20
4	-8	20	74	0.05	250	39	6.5	20
5	-7.5	20	75	0.1	150	40	7	20
6	-7	20	76	0.15	75	51	7	100
7	-6.5	20	77	0.2	40	52	7.25	100
8	-6	20	25	0.25	20	53	7.375	100
9	-5.5	20	26	0.5	20	54	7.5	100
10	-5	20	79	0.525	40	55	7.625	100
11	-4.5	20	80	0.55	75	56	7.75	100
12	-4	20	81	0.6	150	57	7.875	100
13	-3.5	20	78	0.625	300	58	8	100
14	-3	20	82	0.65	150	59	8.125	100
15	-2.5	20	83	0.7	75	60	8.25	100
16	-2	20	84	0.725	40	68	8.375	100
17	-1.75	20	27	0.75	20	61	8.5	100
18	-1.5	20	28	1	20	62	8.75	100
19	-1.25	20	29	1.5	20	63	9	100
20	-1	20	30	2	20	64	9.25	100
21	-0.75	20	31	2.5	20	65	9.5	100
22	-0.5	20	32	3	20	66	9.75	100
23	-0.25	20	33	3.5	20	67	10	100
70	-0.2	40	34	4	20

ID	Center	$k$	ID	Center	$k$	ID	Center	$k$
2	-7	20	59	-1.375	100	30	4	20
3	-6.5	20	60	-1.25	100	31	4.5	20
4	-6	20	61	-1.125	100	49	4.5	100
5	-5.5	20	44	-1	100	50	4.625	100
6	-5	20	45	-0.9375	100	51	4.75	100
7	-4.5	20	46	-0.875	100	52	4.875	100
8	-4	20	47	-0.8125	100	32	5	20
9	-3.5	20	20	-0.75	20	53	5	100
10	-3	20	48	-0.75	100	33	5.5	20
11	-2.75	20	21	-0.5	20	34	5.75	20
12	-2.5	20	22	0	20	35	6	20
13	-2.25	20	23	0.5	20	36	6.5	20
14	-2	20	24	1	20	37	7	20
54	-2	100	25	1.5	20	38	7.5	20
55	-1.875	100	26	2	20	39	8	20
56	-1.75	100	27	2.5	20	40	8.5	20
57	-1.625	100	28	3	20
58	-1.5	100	29	3.5	20

Well Tempered Metadynamics

In addition to REUS, we performed experiments with well-tempered metadynamics to support our REUS results (Figure S11 and Figure S12). We found barrier values within error of our REUS measurements. Because WT MetaD took much longer to converge for such a long RC trajectory, we added soft harmonic wall potentials along the reaction coordinate in both the acylation and deacylation simulations to limit the sampling of only the first transition state (TS1). For the acylation simulations, lower and upper walls were placed at $\mathrm{RC}_{\mathit{acylation}} = (-11.0)$ and 2.2, respectively, with a force constant of 200 kJ mol$^{-1}$. Similarly, for the deacylation simulations, lower and upper walls were applied at $\mathrm{RC}_{\mathit{deacylation}} = (-10.0)$ and 2.0, respectively, using the same force constant. In addition, an upper harmonic wall was applied to the distance between the attacking water oxygen atom and the substrate carbonyl carbon ($d_{\mathrm{Wat:O-MHET:C}}$) at 0.35 nm with a force constant of 2000 kJ mol$^{-1}$ to prevent the catalytic water molecule from diffusing away from the active site during the deacylation simulations. We set the initial Gaussian hill height to 4 kJ mol$^{-1}$ and the bias factor to 20. The bias was deposited every 500 steps. WTMetaD simulations were performed until convergence, which happened within 5 ns for both the acylation and deacylation reactions.

Figure S13 shows the time evolution of the bias potential together with the sampling of all collective variables constituting $\mathrm{RC}_{\mathit{acylation}}$. The accumulated bias plateaued within $\sim$5 ns, indicating convergence of the WTMetaD simulation within the sampled region. Furthermore, in both the acylation and deacylation simulations, the constituent CVs were well sampled throughout the region of interest.

NDPK

Simulation Background

A conserved phosphohistidine intermediate forms within this active site, where a divalent cation (typically $\mathrm{Mg}^{2+}$) coordinates the nucleotide phosphates, balancing their negative charge and organizing the phosphoryl-transfer geometry. As such, accurate modeling of the dephosphorylation half-reaction is sensitive to the local electronic structure, cation coordination, and the local solvent context. The phosphoryl transfer step has proven difficult to model with conventional methods. Cluster models of the active site typically neglect the surrounding protein, ions, and solvent [11].

Simulation Setup and Details

Structurally, NDPK is a compact phosphoryl-transfer enzyme ($\sim150$ amino acids per subunit) that assembles into a homohexamer, with each subunit containing a conserved catalytic histidine. NDPK enables the conversion of ATP to ADP, preceding phosphoryl transfer to an incoming nucleoside diphosphate (NDP). High-resolution crystal structures of transition-state-like complexes such as PDB 1KDN have captured phosphoenzyme and transition-state analog geometries, providing unusually direct structural constraints on the catalytic mechanism. In PDB 1KDN, Dictyostelium NDPK is captured with $\mathrm{Mg}$–ADP and aluminum fluoride, where $\mathrm{AlF}_{3}$ mimics the transferring phosphoryl group [65]. The $\mathrm{AlF}_{3}$ moiety adopts a trigonal-planar arrangement that bridges the leaving-group oxygen of ADP (the $\beta$-phosphate O7) and the N$\delta$ of the catalytic histidine (HIS122), consistent with partial bonding to both the nucleophile and the leaving group in an in-line $\mathrm{S_N2}$-like transition state.

The initial simulation box was prepared using the CHARMM-GUI input generator [105] by placing one NDPK hexamer with ATP docked inside the reactive pocket from the 1KDN [106] crystal structure within a cubic box with an edge length of 7 nm. We employed periodic boundary conditions in all dimensions. The residues were protonated for a pH of 7 except for HIS122 where we enforced a doubly-protonated state as suggested by Hutter and Helms [11]. The system was then solvated in water to a density of 1 g/cm$^{3}$. The total charge of the system was neutralized by adding a potassium ion. To generate the reactant state, consisting of the ATP substrate docked in the NDPK binding pocket, the $\mathrm{AlF}_{3}$ residue in the 1KDN structure was replaced by a $\mathrm{PO}_{3}$ residue and the system was energy minimized using the L-BFGS algorithm [107] until the maximum force was below 0.05 eV/Å to eliminate high-energy overlaps. Equilibration was first performed for 100 ps in the NVT ensemble at 300 K using a Langevin thermostat with a friction constant of 1 ps$^{-1}$. Initial velocities were assigned from a Maxwell–Boltzmann distribution at 300 K. The classical equations of motion were integrated using a 0.5 fs time step. This was followed by 300 ps of equilibration in the NPT ensemble at 300 K and 1 bar using the isotropic Martyna-Tobias-Klein integrator [108] with a temperature coupling time constant of 0.1 ps and a pressure coupling time constant of 1 ps. Final scale runs were conducted in the NVT ensemble at 300 K employing a Langevin integrator with a friction constant of 1 ps$^{-1}$. All MD simulations were conducted on 64 NVIDIA A-100 GPUs, yielding a performance of $\sim$1 ns/day.

In our simulations, the bound ATP adopts hydrogen-bonding contacts to TYR56, ARG92, THR98, and ASN119 in agreement with the 1KDN crystal structure [106] and prior semi-empirical results [11], and the $\mathrm{Mg}^{2+}$ ion spontaneously adopts octahedral hexacoordination with the three phosphoryl groups of ATP, two water molecules, and GLU58 (Figure 4B). This geometry emerges without ion-specific empirical corrections. In contrast, MM approaches often fail to reproduce the $\mathrm{Mg}^{2+}$ solvation free energy and first-shell geometry [109, 110], and QM/MM simulations require a priori specification of the residues, solvent molecules, and ions participating in the QM region.

Enhanced Sampling Calculations

Enhanced sampling of the phosphoryl transfer reaction was conducted to study the reaction mechanism and estimate characteristic transition times. To understand the thermodynamics and the mechanism of the phosphotransfer reaction, we used two collective variables (CVs) which track nucleophilic attack of HIS and the proton transfer from the doubly-protonated HIS to the phosphoryl group. Both CVs were constructed using coordination numbers [111]. We define the coordination number $S$ between the atom groups $A$ and $B$ as,

where,

In the above definition, $d_0$ is the minimum inter-particle distance, $r_0$ is the midpoint of the switching function curve, and $n$ and $m$ are non-negative numbers which control the steepness of the switching function curve. We constructed the first collective variable, CV$_1$, to track the phosphoryl transfer process as,

where $S_{\mathrm{P-N}}$ represents the coordination number of $\gamma$-P of the phosphoryl group with the acceptor N$\delta$ on HIS and $S_{\mathrm{P-O}}$ refers to the coordination number of $\gamma$-P with all the oxygen atoms on the $\beta$-phosphate groups of ADP. Similarly, the second collective variable CV$_2$ tracks the proton hopping process from HIS to the transferring $\gamma$-phosphoryl group,

where $S_{\mathrm{O-H}}$ corresponds to the coordination number of the proton on N$\delta$ of HIS with all oxygen atoms on the transferring $\gamma$-phosphoryl group, and $S_{\mathrm{N-H}}$ represents the coordination number of the proton with the N$\delta$ atom on HIS. Note that the inclusion of all oxygen atoms in the $\beta$-phosphate groups of ADP in $S_{\mathrm{P-O}}$ and all oxygen atoms in the transferring $\gamma$-phosphoryl group in $S_{\mathrm{O-H}}$ automatically includes all degenerate states.

The CVs defined above were used in a 2D WTMetaD simulation to study the free-energy landscape of the phosphotransfer reaction and extract mechanistic insight. During the WTMetaD simulations, a flat-bottom restraint was applied to the sum of CV$_1$ and CV$_2$ to prevent the proton from migrating away from the reaction site via water-mediated hops. A one-sided flat-bottom restraint was imposed between the center of mass of ATP and HIS122 to prevent the ATP substrate from leaving the active site. In practice, ATP remained docked throughout all calculations and the restraint never exerted a restoring force. A convergence assessment of the WTMetaD calculations is presented in Figure S15.

We further performed MD simulations with umbrella sampling along the phosphoryl transfer pathway using CV$_1$ and along the proton transfer pathway using CV$_2$ to improve the estimate of the location and height of the barrier along each path. The results of the MD simulations with umbrella sampling were reweighted using the BayesWHAM algorithm [92] to calculate the potential of mean force (PMF) curve along CV$_1$ and CV$_2$ and their relative free energy. For the umbrella sampling runs on CV$_1$, we use harmonic restraints with a force constant of 150.0 kJ mol$^{-1}$ to keep CV$_2$ at 0.75. Umbrellas were placed along the phosphoryl-transfer coordinate and their force constants were tuned such that there was sufficient overlap between neighboring windows and the drift from each window was within acceptable range. Table S7 shows the centers and force constants of the final umbrella setup. The samples collected from each window were then reweighted using the BayesWHAM algorithm [92] to get the relative free energy of the windows and calculate the final PMF. A convergence assessment of these US calculations is presented in Figure S16. We repeated this process with umbrellas along $\mathrm{CV}_2$, using harmonic restraints with a force constant of $150.0\ \mathrm{kJ}\ \mathrm{mol}^{-1}$ to keep $\mathrm{CV}_1$ at (-0.65), to get an estimate of the proton hopping barrier along the H transfer coordinate. Table S8 shows the centers and force constants of the final umbrella setup. Similar to before, the samples collected from each window were then reweighted using the BayesWHAM algorithm [92] to get the relative free energy of the windows and calculate the final PMF. A convergence assessment of these US calculations is presented in Figure S17.

Index	Center	Force const. ($\mathrm{kJ}\ \mathrm{mol}^{-1}$)
1	-0.75	1000.00
2	-0.60	1000.00
3	-0.45	1000.00
4	-0.35	2000.00
5	-0.35	2750.00
6	-0.30	1000.00
7	-0.25	2500.00
8	-0.15	1000.00
9	0.00	1000.00
10	0.00	2000.00
11	0.15	1000.00
12	0.25	2000.00
13	0.30	1000.00
14	0.45	1000.00
15	0.60	1000.00

Index	Center	Force const. ($\mathrm{kJ}\ \mathrm{mol}^{-1}$)
1	-0.750	1000.00
2	-0.500	1000.00
3	-0.425	2000.00
4	-0.250	1000.00
5	-0.150	2000.00
6	0.000	1000.00
7	0.150	3000.00
8	0.250	3000.00
9	0.375	3000.00
10	0.500	3000.00

Marginalization of the Reweighted WTMetaD Samples

To compare the FES obtained from the 2D WTMetaD run with the umbrella sampling results, we marginalized the 2D surface along CV$_1$ and CV$_2$. Figure S18 shows the marginalized PMFs superimposed on the umbrella sampling results. The location and height of the barriers along both CVs are in good agreement. We observe that the product of the H transfer reaction along CV$_2$ is at a higher free energy from the umbrella sampling results. We attribute this to the difference between how the product region is constrained between the two PMFs: in the umbrella sampling runs, a restraint on the P-O bond prevents the formation of the final product, whereas in the marginalized results the samples collected for CV$_1 > (-0.3)$ are excluded to isolate the destabilizing effect of the transition state as well as the final product along CV$_1$. Note that the PMF from US shown here is calculated without block-averaging.

Cumulative Distribution of Reweighted OPESf Simulations for NDPK

We applied OPESf to the rate limiting phosphoryl transfer step to calculate rate constants and effective activation energies using Eyring transition-state theory (Figure 4E, Section A.4) [49, 72]. Reweighted first-passage times from 19 independent OPESf trajectories follow Poisson statistics, yielding a predicted rate constant of $k_{\mathrm{OPESf}} = (0.043 \pm0.019)\ \mathrm{s}^{-1}$ with a corresponding apparent effective activation free energy of $E_{a}^{\mathrm{OPESf}} = (81.32 \pm1.13)\ \mathrm{kJ}\ \mathrm{mol}^{-1}$ (Figure S19).

Water Structure and Dynamics

We evaluated eSEN-omol’s capabilities in predicting the temperature dependence of the structural, thermodynamic, and kinetic properties of bulk water. An initial box of water containing 216 molecules was generated using GROMACS [112] and geometrically relaxed with the SPC/E water model [113]. Molecular dynamics simulations with eSEN-omol were conducted using the integrators implemented in ASE [81]. All calculations were performed at neutral charge and in the singlet spin state.

We first computed radial distribution functions (RDFs) as a function of temperature to assess the structural behavior of eSEN-omol for liquid water. These RDFs were calculated from equilibrated trajectories of 250 ps in the NVT ensemble using the Langevin integrator with a low friction coefficient of $0.001\ \mathrm{fs}^{-1}$ and a time step of $0.5\ \mathrm{fs}$, employing the eSEN-omol model. Twelve independent simulations were performed at temperatures ranging from 248 K to 368 K. The resulting RDFs were compared against experimental X-ray diffraction measurements, for which we observe very good agreement [114, 115] (Figure S20).

We next used these simulation trajectories to evaluate the self-diffusion coefficient, $D_\mathrm{s}$, from the mean-squared displacement (MSD) using the Einstein relation, $D_\mathrm{s} = \lim_{t \to\infty} \frac{\langle\lvert\mathbf{r}(t)-\mathbf{r}(0) \rvert^2 \rangle}{6t}$,

where the MSD was computed independently for each oxygen atom and averaged over all 216 atoms, with a linear fit applied in the diffusive regime (lag times 2–50 ps). It is well known that simulations employing periodic boundary conditions systematically underestimate the self-diffusion coefficient because a molecule couples hydrodynamically to its own periodic images through the long-ranged Oseen flow field [116]. We applied the Yeh–Hummer finite-size correction, $D_\infty= D_\mathrm{PBC} + \frac{k_\mathrm{B}T\xi}{6\pi\eta L}$,

where $\xi= 2.837297$ is a dimensionless geometric constant for a cubic periodic box, $L$ is the box side length, and $\eta$ is the shear viscosity of water. Reference viscosity values were taken from the IAPWS international formulation for $T \ge273$ K [117] and from microfluidic Brownian-motion measurements in the supercooled regime for $T < 273$ K [118]. The corrected diffusivities predicted by eSEN-omol are in excellent agreement with MB-pol [119] and experimental measurements [120, 121, 122, 123] (Figure S21).

Finally, we evaluated the density of liquid water at 300K and 1 atm pressure. We generated boxes of water with TIP3P geometry using Open-MM (1k and 10k waters) and ran isotropic Martyna–Tobias–Klein [108] with a temperature damping of 100 fs and a pressure damping of 1000 fs for 200 ps each with 0.5 fs time steps. eSEN-omol, like all other OMol-trained models, is known to have an elevated water density due to the neglect of many-body dispersion in the VV10 functional used as part of the $\omega$B97M-V functional in O​​Mol25 [32, 124, 125]. We show that we can correct this problem by fine-tuning eSEN-omol with a similar functional where 3-body dispersion is included. The OMol-4M subset [32] was recomputed at $\omega$B97M-D4, using the D4 parameters of Friede et al. [126] (this is known as “$\omega$B97M-D4rev” in ORCA[127] and includes the 3-body ATM correction). Results are shown in Figure S22.

In Figure S23, we show that using the D4-corrected model, we can closely recover the barrier height compared to the uncorrected model. We also show that we can run all our simulations close to 1 g/cm$^3$ densities using NVT alone with no NPT using the uncorrected eSEN-omol model to produce nearly identical results. This result is not very surprising: the neglect of 3-body dispersion introduces an extremely systematic shift in the absolute interaction energies of molecules whose effect becomes manifest mainly when these small quantities accumulate at the condensed phase scale; relative energies, such as barrier heights, are expected to be only barely affected.

Folding of Polyalanine in Vacuum

The polyalanine peptide with a length of 15 residues (ACE-Ala15-NME) is known to form a stable helix structure at 300 K in vacuum [128, 129]. We used this peptide to assess eSEN-omol in forming the expected stable helix structure from an extended configuration. We constructed the initial configuration of ACE-Ala15-NME in an extended configuration using AmberTools v24.8 [130]. We then minimized the energy of this structure using the LBFGS algorithm in ASE [81], setting a convergence criterion for the force on all individual atoms to be less than 0.01 eV/Å within a maximum of 10,000 steps. Then, we conducted Langevin dynamics for 1 ns with a friction coefficient of 0.01 fs$^{-1}$ and a time step of 1 fs in a vacuum at 300 K using the eSEN-omol model s1p1 calculator in ASE [81]. The net charge and spin of the system were set to 1 and 0, respectively. The changes in the secondary structure of ACE-Ala15-NME during the 1 ns simulation, as defined by the Dictionary of Secondary Structure in Proteins (DSSP) [131], are shown in Figure S24. Within 200 ps, the extended coil configuration transitions into a bend configuration and then into a $3_{10}$ helix structure. The predominant $3_{10}$ helix structure lasts up to 0.7 ns before $\sim$75% of residues in Ala15 convert to the stable $\alpha$-helix structure. A few of the residues form a turn and retain the coil configuration during the last 0.2 ns of the simulation. The observed variations between the different types of secondary structures with a stable $\alpha$-helix structure at the end of the 1 ns simulation are in good agreement with [128] and [129].