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 remains 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]. Computational expense can be mitigated with semi-empirical QM methods but at a cost in accuracy [16, 17]. When solvent plays a role in the reaction, the manually chosen number QM water molecules adds another parameter that can strongly influence results [18, 19]. QM/MM still remains computationally expensive, with simulations typically limited to tens of picoseconds, [19, 20] 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 AI$^2$BMD, which achieves ab initio protein energies and forces using a fragmentation strategy and treats solvent using a polarizable model [21]. 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 [22, 23]. 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 [24, 25, 26], but the challenges associated with ML region specification, MM force field choice, ML/MM boundary treatment, and solvent modeling all remain.

Recent reports have shown that models trained on Open Molecules 2025 dataset [28] match DFT performance in numerous chemical and biological applications [25, 29, 30, 31].

We show here that eSEN-omol¹, the base model of UMA [32], 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 [28]. We leverage established enhanced sampling protocols for rare-event simulation to promote transition state crossing and accurately estimate barrier heights in enzymatic reactions [33]. 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 A.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 demonstrating that it is able to reproduce sampled cluster DFT single-point energies with high-fidelity, correctly predict long-range electrostatic behavior in scenarios relevant to proteins in solvent environments, and keep proteins correctly folded over many nanoseconds of simulation. Figure 1B-E shows a collection of these validations with supporting details provided in Appendix A.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 and demonstrate the accuracy and transferability of our method using the same eSEN-omol potential with no additional parameterization.

Our results demonstrate a unified approach to accurately modeling free energy barriers and reaction rates in complex enzyme-catalyzed reactions, establishing a practical 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 [34, 35]. 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 [36]). Experimentally, the enzyme accelerates the reaction by $10^6$ to $10^{-7}$-fold relative to aqueous solution [37], 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 [35, 38, 39].

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) [40, 41] along a previously reported collective variable (CV) to drive the Claisen rearrangement (Figure 2B, Appendix B) [19, 34]. The resulting free-energy surface (FES) in water gives a barrier of $\Delta G_{\mathrm{water}}^{\ddagger} = (119.81 \pm6.24)$ kJ mol$^{-1}$ (Figure 2) in good agreement with experiment ($102.55$ kJ mol$^{-1}$ [42]) and previous computational studies

(102.09 $\pm$ 2.51 kJ mol$^{-1}$ [18]) (fig: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 \pm2.72)$ kJ mol$^{-1}$, in reasonable agreement with the experimentally-reported value of 64.43 kJ mol$^{-1}$ [37] (fig:2C,fig:2D), and corresponding to a $\sim45$ kJ mol$^{-1}$ reduction compared to the barrier in water. The $\sim10$ kJ mol$^{-1}$ discrepancy with respect to experiment is comparable to the errors in the DFT training energies underpinning eSEN-omol [43]. 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 $\sim6.5$ kJ mol$^{-1}$.

We additionally computed rate constants using OPES flooding (OPESf) [45]. In water, we obtained a rate constant of $k_{\mathrm{water}} = (9.01 \pm4.53) \times10^{-9}\ \mathrm{s}^{-1}$ corresponding effective activation free energy of $E_{\mathrm{a}}^{\mathrm{water}} = (119.69 \pm1.28)\ \mathrm{kJ\ mol}^{-1}$ (Figure S7). In the enzyme, the calculated rate constant increases to $k_{\mathrm{enzyme}} = (0.26 \pm0.12)\ \mathrm{s}^{-1}$ implying an effective activation free energy of $E_{\mathrm{a}}^{\mathrm{enzyme}} = (76.77 \pm1.16)\ \mathrm{kJ\ mol}^{-1}$. The calculated kinetic activation energies are in excellent agreement with the free energy barriers within $\sim1.5\ \mathrm{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 [37]. These results constitute a marked improvement on replicating experimental barriers over previous ML/MM [44] and QM/MM [19] studies at comparable (high) electronic-structure levels of theory (fig:2D). We also note that the QM/MM study found that QM region selection had a profound effect on the computed barrier [19], a complication which eSEN-omol avoids.

Ester-Bond Hydrolysis in PETase

The PETase enzyme from Ideonella sakaiensis is a promising biocatalyst for depolymerization of poly(ethylene terephthalate) (PET) [46]. PETase employs the canonical catalytic SER-HIS-ASP triad for PET hydrolysis [47, 48] (fig:3A, fig:3B). Structural and biochemical studies suggest that PET hydrolysis proceeds through sequential acylation and deacylation steps [46]. 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 [49, 50, 51, 52, 53, 54], but differences in DFT settings [50, 53], selection of the QM region [53, 55], length of the polymer chains [51], and choice of CVs [49, 53] has led to debate on whether the acylation [20, 50] or deacylation [49, 56] step is rate limiting and whether each reaction proceeds in a single step [49] or through a two-step tetrahedral intermediate [51].

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

As illustrated in fig:3C, we predict an effective free energy barrier of $(79.03 \pm4.2)\ \mathrm{kJ\ mol}^{-1}$, derived from the acylation step, to be rate limiting, in good agreement with prior QM/MM studies [20, 49, 50, 51] and experimental measurements for the BHET dimer [46]. We assume, as have others [20, 49, 50, 51], that the acylation and deacylation step 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 heights of $(75.19 \pm0.4)\ \mathrm{kJ\ mol}^{-1}$ [46, 59] and that predicted by Burgin et al. of $76.5\ \mathrm{kJ\ mol}^{-1}$ [49]. 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 [20, 50]. The highest free energy along the acylation reaction at transition state TS2 is $(79.03 \pm4.2)\ \mathrm{kJ\ mol}^{-1}$ (fig:3D, Table S5), whereas the highest barrier along the deacylation reaction at transition state TS3 is slightly lower at $(60.72 \pm3.0)\ \mathrm{kJ\ mol}^{-1}$ (fig:3E, Table S5). Second, despite using the identical CV and similar umbrella sampling methods, we predict both acylation and deacylation to proceed by a two-step mechanisms with well-defined intermediates as opposed to a single concerted step. The acylation intermediate I1 lies $\sim$19 $\mathrm{kJ\ mol}^{-1}$ lower in free energy than TS1 and $\sim$25 $\mathrm{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 $\mathrm{kJ\ mol}^{-1}$ lower in free energy than TS3 and $\sim$19 $\mathrm{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 tetrahedral geometry centered on the PET ester carbonyl carbon, consistent with the simulations of Guo et al. [51], which employed a higher level of DFT theory (M06-2X/6-31G) than Burgin et al. (DFTB3) [49] 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 [51, 53]. Taken together, our results suggest that eSEN-omol, which is trained at a $\omega$B97M-V/def2-TZVPD level of theory, can accurately predict barrier heights and reaction mechanisms at

Reaction	Feature	Mean RC	$\Delta G$ (kcal $\mathrm{mol}^{-1}$)	95% CI (kcal $\mathrm{mol}^{-1}$)	$\Delta G$ (kJ $\mathrm{mol}^{-1}$)	95% CI (kJ $\mathrm{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]

∼1000× speedups [19, 51] relative to equivalent QM/MM methods.

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 [60, 61].

It operates through a ping–pong mechanism in which the active-site HIS122 captures the $\gamma$-phosphoryl group of a nucleoside triphosphates (NTP) such as ATP, and donates it to an incoming nucleoside diphosphate (NDP) [62, 63]. It has been shown that key energetics of QM/MM approaches often fail to converge until the QM region is expanded to hundreds of atoms [8, 64, 65, 66], making simulation exceedingly costly and limited to short time scales. Therefore, we sought to 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 [60]).

We calculated the FES of the phosphoryl transfer reaction using well-tempered metadynamics (WTMetaD) [68] driven in two collective variables: CV$_1$ tracks nucleophilic attack of HIS122 $N\delta$ on the $\gamma$-phosphorus, and CV$_2$ tracks proton transfer from $N\delta$ of HIS122 to the $\gamma$-phosphoryl oxygen (Appendix D). 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 S$_\mathrm{N}$2-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 ∼7.5 kJ mol$^{-1}$. The transition state (VII) of the concerted proton/phosphoryl transfer mechanism (I→VII→V) stands ∼100 kJ mol$^{-1}$ above the reactant, which is substantially higher than the corresponding value of ∼70 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 OPESf rate calculation.

As an independent verification of the free energy barriers calculated from 2D WTMetaD, we conducted 1D umbrella sampling (US) [69] of the two elementary steps in the proposed two-step pathway. US calculations in CV$_2$ between states I and III exposes 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 CV$_1$ between states III and V resolves 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 S21). 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 OPESf to the rate limiting phosphoryl transfer step (Figure 4E, Section A.4) [45, 70] and obtained a predicted rate constant of $k_\mathrm{OPESf} = (0.043 \pm0.019)$ s$^{-1}$ with a corresponding apparent effective activation free energy of $E_a^\mathrm{OPESf} = (81.32 \pm1.13)$ kJ mol$^{-1}$ (Figure S22). The effective activation energy agrees with the free energy barrier computed by US to within ∼10 kJ mol$^{-1}$. In Figure 4G we compare our US (73 kJ mol$^{-1}$) and OPESf (80.3 kJ mol$^{-1}$) barriers to the value implied by the experimental $k_\mathrm{cat}$ (55.6 kJ mol$^{-1}$ [67]). 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 ∼17 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

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 at a $\omega$B97M-V/def2-TZVPD level of theory, enables higher accuracy than conventional QM/MM calculations which often employ the B3LYP or PBE functionals that are known to underestimate reaction barriers, and even relatively more accurate functionals like the rarely employed M06-2x and $\omega$B97X-3c [71, 72]. Additionally, eSEN-omol offers $\sim1000\times$ speedups over QM/MM, 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 [58], WTMetaD [68], OPES [41], and OPESf [45] – 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 $\mathrm{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 [28] 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 to not only understand complex biomolecular processes, but 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. 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

References

[1] [2] [3] [4] [5] [6] [7] [8]

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. 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], are limited to inter-node only parallelism [87], or require explicit compatibility with the LAMMPS [88, 89] interface which is often not compatible with modern MLIPs. To the best of our knowledge, no existing MLIP can run MD at speeds of nanoseconds/day on large solvated protein systems with $O(100\mathrm{k})$ atoms. 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

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. Section 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-position active-site cluster was constructed from each configuration. Each 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 [90, 91], 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 (fig:1C, Figure S2, Figure S8, and Figure S17). 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 [92].

Evidence for Correct Long Range Behavior

Highly charged peptide. A 2 ns simulation of the highly charged 1FME peptide [93] comprising 14 charged residues (out of 28 total) and numerous long-range salt bridges produces a close match to experimental NMR structures (fig: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 (fig:1D). Replicating the calculations from Ref. [27], we perform NVT simulations of a water box with 512 atoms over 2 ns and show the long range dipole correlations function ($\langle\tilde{m}^{*}(k)\tilde{m}(k)\rangle$, where $\tilde{m}$ is the Fourier transform of the molecular dipole density, agree 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 [27], we use the same box size with 512 atoms and 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.

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

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 an NaCl ion pair closely agrees with force fields with explicit long range electrostatics terms (Figure 1). 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 1).

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 2). Moreover, eSEN-omol’s PMF, fitting the coefficients of an $A/r + C$ functional form over 8–18 Å 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 [94], 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 Appendix fig:E and fig:F, 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 Columbic 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 [95]. 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 [95, 96, 97, 98, 99]. 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 [99, 100, 101].

Uncertainties in the OPESf transition times and rate constants were estimated following the procedure described in Ref. [70]. 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. [102], with the protein conformation based on PDB 3ZO8 [36]. 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 [32]. 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 [19], the Claisen rearrangement of chorismate to prephenate was described using a single CV defined as 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 CV = (−0.05) nm and an upper wall at 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 (fig:S3a–fig:S3c), satisfying the multiple-crossing criterion for sampling convergence [41, 103]. The number of deposited OPES kernels plateaus at ∼111 (fig:S3d), 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; fig:S3e–fig:S3f) 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 (fig:S3g).

Going beyond equilibrium free energies, we next computed rate constants using OPES flooding (OPESf) [95]. 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) [45, 70]. In water, we computed a rate constant of $k_{\mathrm{water}}=(9.01 \pm4.53)\times10^{-9}$ 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)$ 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 [37].

Structural Model Variation

While the prepared enzyme-substrate structures used in Refs. [19, 44] were not provided, we note that when we attempted to dock chorismate into the PDB: 2CHT [104] 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). In contrast, the higher-resolution PDB: 3ZO8 [36] structure resulted in a more stably docked substrate, featuring a two-hydrogen-bonded substrate conformation and a substantially lower activation barrier. 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.

PETase

Simulation Setup and Details

Initial structures for both the acylation and deacylation simulations were obtained from Burgin et al. [49]. The deacylation simulations were initiated from the same acyl-enzyme intermediate (AEI) structure used in Burgin et al. [49]. In all the PETase simulations, we employed ASE [81] together with eSEN-omol [32]. 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}/\text{\AA}$. 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 Burgin et al. [49] for both the acylation and deacylation reactions. The acylation reaction coordinate ($\mathrm{RC}_{\mathrm{acylation}}$) is a linear combination of three collective variables (CV1a, CV2a, CV3a). CV1a and 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. 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 (CV1d, CV2d, CV3d). CV1d, CV2d, and 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

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) [58]. In this variant of umbrella sampling (US) [69], 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 criteria 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 [40] 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 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 S3 and Table S4, 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 [94, 106]. We detail the key features on the free energy surfaces in Table S5. The total stable sampling time was $\sim400$ ps per window, total about 50 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 validate our results (Figure S11 and Figure S12). We found barrier values within error of our REUS measurements. Because WT MetaD was harder 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. For the acylation simulations, lower and upper walls were placed at $RC_{\mathrm{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 $RC_{\mathrm{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 $RC_{\mathrm{acylation}}$. The accumulated bias plateaued within $\sim$5 ns, indicating convergence of the WTMetaD simulation. Figure S14 shows the corresponding convergence analysis for the deacylation simulations of both the wild-type and S214H mutant systems (Section C.5). In both cases, the accumulated bias plateaued within $\sim$5 ns, indicating convergence of the WTMetaD simulations. Furthermore, in both the acylation and deacylation simulations, the constituent CVs were well sampled throughout the region of interest.

Altered Aromatic Interactions in the S214H Mutant During the Deacylation Reaction

We analyzed the stacking distance and stacking angle between TRP185 and the aromatic ring of MHET (mono(2-hydroxyethyl) terephthalate) together with the TRP185 side-chain torsion angle ($C - C_{\alpha} - C_{\beta} - C_{\gamma}$). By analyzing simulations of a S214H mutant, we found that the mutation altered the aromatic interactions of TRP185 with MHET (Figure S15), which we hypothesize to be responsible for the increased stability of the AEI, resulting in higher deacylation activation barrier and the experimentally-observed reduced catalytic activity [107]. Our analysis (Figure S15) indicates that the S214H mutant alters the conformational dynamics of TRP185 and potentially strengthens its stacking interactions with MHET, which forms the AEI with SER160. We therefore hypothesize that the altered aromatic interactions in the S214H mutant, reflected by the shorter TRP185–MHET distance suggests that TRP185 adopts conformations that favor $\pi$-stacking with the aromatic ring of MHET over the imidazole ring of HIS214. This contrasts with the WT, where SER185 lacks an aromatic side chain and therefore does not introduce such a competition for aromatic interactions with TRP185. We propose that these altered aromatic interactions are responsible for the increased stability of the AEI, resulting in higher deacylation activation barrier and are consistent with the experimentally-observed reduction in catalytic activity [107].

FES Along Simpler CVs

To test the robustness of our results to the choice of reaction coordinate, we projected the configurations sampled during the WTMetaD simulations onto the CVs employed by Guo et al. [51] as an alternative, simpler CVs with which to characterize the deacylation reaction (Figure S16). Under these CVs, we predict an acylation free energy barrier of ($71.5 \pm5.4$) kJ mol$^{-1}$ and a deacylation free energy barrier of ($66.4 \pm2.4$) kJ mol$^{-1}$, which are in reasonable agreement of the values of ($59.4 \pm0.8$) kJ mol$^{-1}$ and ($68.2 \pm0.4$) kJ mol$^{-1}$ reported by Guo et al. using QM/MM calculations in which the QM region is treated at a M06-2X/6-31G(d) level of theory [51].

NDPK

Simulation Background

A conserved phosphohistidine intermediate forms within this active site, where 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 Mg–ADP and aluminum fluoride, where $\mathrm{AlF}_{3}$ mimics the transferring phosphoryl group [60]. 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 nucleophile and the leaving group in an in-line $\mathrm{S_N2}$-like transition state.

The initial simulation box was prepared using CHARMM-GUI input generator [108] by placing one unit of NDPK hexamer with ATP docked inside the reactive pocket from 1KDN [109] crystal structure within a cubic box with an edge length of 7 nm edge. 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/$\mathrm{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 1KDN structure was replaced by a $\mathrm{PO}_{3}$ residue and the system was energy minimized using the L-BFGS algorithm [110] 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 $\mathrm{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 [111] 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 $\mathrm{ps}^{-1}$. All MD simulations were conducted on 64 NVIDIA A-100 GPUs, yielding a performance of $\sim1$ ns/day.

In our simulations, the bound ATP adopts hydrogen-bonding contacts to TYR56, ARG92, THR98, and ASN119 in agreement with 1KDN crystal structure[109] 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 (fig: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 [112, 113] 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 [114]. 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 atom in the $\beta$-phosphate groups of ADP in $S_{\mathrm{P-O}}$ and all oxygen atom 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 restraints 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 S18.

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 [94] 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 S6 shows the centers and force constants of the final umbrella setup. The samples collected from each window were then reweighted using the BayesWHAM algorithm [94] 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 S19. We repeated this process with umbrellas along $\mathrm{CV}_2$, using harmonic restraints with a force constant of $150.0\ \mathrm{kJ\,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 S7 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 [94] 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 S20.

Index	Center	Force const. (kJ 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. (kJ 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 S21 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 the 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 collect 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) [45, 70]. 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_{\mathrm{a}}^{\mathrm{OPESf}} = (81.32 \pm1.13)\ \mathrm{kJ}\ \mathrm{mol}^{-1}$ (Figure S22).

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 [115] and geometrically relaxed with the SPC/E water model [116]. 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 [117, 118] (Figure S23).

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{\left\lvert\mathbf{r}(t)-\mathbf{r}(0)\right\rvert^{2}}{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 [119]. 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 [120] and from microfluidic Brownian-motion measurements in the supercooled regime for $T < 273$ K [121]. The corrected diffusivities predicted by eSEN-omol are in excellent agreement with MB-pol [122] and experimental measurements [123, 124, 125, 126] (Figure S24).

Finally, we evaluated the density of liquid water at 300 K 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 [111] 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, are known to have an elevated water density due to the neglect of many-body dispersion in the VV10 functional used as part of the ωB97M-V functional in OMol25 [28, 127, 128]. 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 [28] was recomputed at ωB97M-D4, using the D4 parameters of Friede et al. [129] (this is known as “ωB97M-D4rev” in ORCA [130] and includes the 3-body ATM correction). Results are shown in Figure S25 Figure S25.

In Figure S26 Figure S26, 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 who effect becomes manifest mainly when these small quantities accumulate at the condensed phase scale; relative energies, such as in 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 [131, 132]. 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 [133]. 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/$\AA$ 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 simulations, as defined by the Dictionary of Secondary Structure in Proteins (DSSP) [134] are shown in Figure S27. 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 $\sim75\%$ of the fraction 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 Kovács et al. [131] and Kabylda et al. [132].