Application of neural network potentials to modelling transition states

Abstract

Transition state modelling remains a challenge in computational chemistry, often requiring chemical intuition and expensive, iterative recalculations. This work presents a more efficient approach using umbrella sampling to explore free energy surface and more importantly, the conformational space around transition states, reducing the effort needed for structure identification. By employing a machine learning potential, ANI-2x, [C. Devereux et al., J. Chem. Theory Comput., 2020, 16, 4192–4202] to drive the sampling, we demonstrate enhanced FES exploration and efficiency compared to traditional DFT methods. The approach is applied to two different reactions: amide formation via a thioester intermediate and disulphide bridge formation. It was found that ANI-2x performs poorly at the prediction of high energy structures yet provides rapid, thorough sampling of reaction pathways making it useful for informing further calculations at higher levels of theory.

---

The exploration of chemical reaction pathways has long been a cornerstone of advancing our understanding of chemistry. This pursuit has given rise to techniques such as Density Functional Theory (DFT) and quantum chemistry, which have become essential tools for underpinning the exploration of structure and reactivity. More recently, machine learning potentials (MLPs) have made large strides in the efficient and accurate calculation of atomic energies and forces. MLPs represent a middle point between DFT and molecular mechanics (MM), combining the accuracy of DFT and other quantum methods with the speed of MM. MLPs are trained to learn the shape of the potential energy surface from reference data, e.g., derived from DFT calculations.

In recent years, several different MLP architectures have been published, capitalizing on different machine learning techniques such as neural networks. Recently, ANI-2x extended support to 7 atom types (H, C, N, O, F, S, Cl) and boasts accuracy comparable to DFT methods but 10⁶ times faster.^1,2

Other MLPs have since been developed with varying architectures such as MACE and AIMNet and have been applied to a wide range of applications.^3,4 These applications range from the calculation of chemical properties such as pK_a values⁵ to running MD simulations that extend to hundreds of nanoseconds for systems with a large number of atoms,⁶ as well as studying bulk systems of substances such as water, metals, and perovskites.⁷

A particularly prevalent area of application for integrated ML techniques is the study of chemical processes and reactions. Reactivity within explicit solvent,⁸ catalytic reactions on surfaces,⁹ nucleation structures,¹⁰ and pericyclic reactions¹¹ have all benefitted from the acceleration that MLPs offer over traditional quantum chemical methods, whilst retaining some of their accuracy. Accurate modelling of such systems with traditional static quantum chemistry approaches (single minimum and transition state) offers limited insight and often does not apply at higher temperatures where increasing numbers of conformations can be accessed. By instead pursuing extensive conformational sampling through techniques like umbrella sampling (US),¹² we can provide a more holistic view of molecular systems.

In this work, the ANI-2x MLP is employed to drive an Umbrella Sampling based workflow, studying two different reactions: the formation of an amide (reaction 1) and the formation of a disulphide bridge (reaction 2) as shown in Scheme 1. These reactions were chosen as they represent key chemical processes in large biochemical systems. Such large-scale systems are not tractable using pure DFT approaches and often require simplifications or hybrid methods. Therefore, the ability of ANI-2x to accurately model these processes within complex systems would be highly informative. Further, we explore ANI-2x's ability to (a) reproduce DFT energies for given stationary point structures, and (b) explore the free energy surface for a reaction through enhanced sampling methods.

Scheme 1 Reactions studied: (1) peptide bond formation through thioester intermediate with H₂S loss. (2) Disulphide bridge formation via hydrogen peroxide splitting.

Firstly, the reaction mechanisms for reactions 1A, 1B and 2, were explored with DFT and the ANI-2x potential. DFT calculations were performed at the ωB97X/6-31G* level of theory^13,14 in ORCA V6.0.1,^15,16 to ensure consistency with the ANI-2x dataset.¹⁷ Note that this level of theory is suitable for reaction pathway exploration rather than definitive energetic benchmarking. While ωB97X does not accurately capture dispersion, this was found not to be a significant contribution to the reaction study, ESI† (page S5). The ANI-2x model was implemented in Python through the TorchANI module¹⁸ and was used through an atomic simulation environment (ASE) calculator for single point calculations.¹⁹ As shown in Fig. 1, the formation of the intermediate thioester species would be unlikely to proceed at room temperature, consistent with previous reports requiring catalysis.²⁰ Additionally, it should be noted that the C-terminus sulphur species in native chemical ligation is usually reported as a C-terminus thioester²¹ rather than the capped thiol species we use here to specifically limit our system size and complexity and as such this may influence the resultant energy barriers shown.

Fig. 1 Relative energies (kcal mol⁻¹) for DFT (B97X/6-31G*, black) and ANI-2x (red) relative to reactants.

ANI-2x shows varying accuracy when predicting high-energy structures. For transition states in reactions 1A and 1B, ANI-2x exhibits differences of 10.9 kcal mol⁻¹ and 3.9 kcal mol⁻¹ respectively when compared to DFT calculations, values which clearly exceed accepted chemical accuracy of 1 kcal mol⁻¹. These findings align with previous results that found that ANI-2x overestimates proton transferred structures and fails to generalise well to high energy isomers.¹⁰ Contrastingly, ANI-2x demonstrates reasonable performance for the equilibrium species in reaction 1, with differences of only 1.5 and 1.2 kcal mol⁻¹ for intermediate and product structures respectively. However, ANI-2x's accuracy varies in certain chemical environments. In reaction 2, the 4.9 kcal mol⁻¹ difference for the product species exceeds chemical accuracy (ESI,† Fig. S2).

Reaction 2 was explored in a similar manner to reaction 1 via DFT. The transition state barrier energy was found to be 80.1 kcal mol⁻¹ at the ωB97x/6-31G* level of theory. This computational method was selected to maintain consistency with the ANI-2x dataset, though we also explored other theoretical approaches as detailed in the ESI† (pages S2 and S3). Previous studies have observed variations up to 14.1 kcal mol⁻¹ between different functionals using the 6-311G basis set for similar reactions.²² Our evaluation across different levels of theory suggests that reaction 2 exhibits self-interaction error (SIE), which can be partially addressed by increasing HF exact exchange. For this reaction, ANI-2x calculates a product species relative energy of −71.6 kcal mol⁻¹ compared to DFT's −66.7 kcal mol⁻¹ (ESI,† Fig. S2). Additionally, ANI-2x predicts a forward energy barrier of 17.2 kcal mol⁻¹, which is a 62.9 kcal mol⁻¹ difference from the DFT value.

This difference likely reflects ANI-2x's training dataset, which lacks systems with SIE or transition states involving concurrent bond breaking/forming. Notably, ANI-2x predicts equilibrium species energies more accurately than transition state energies.

So far, we have detailed the ability of ANI-2x in predicting structures with reference to potential energies of stationary point structures. However, it is possible to optimise structures in a similar approach to DFT (as shown in the ESI,† Section S3) although this still limits the exploration of reaction space. Further, an intrinsic part of computational chemistry relies upon the thermodynamic properties of systems. ANI-2x is trained on electronic energies and forces, however, we deploy ANI-2x here for the calculation of forces in US simulations, predicting free energies and monitoring systems via enhanced sampling along reaction coordinates (RCs). Details of RC_1A, RC_1B and RC_2, the RC for each respective reaction can be found in Table 1 with atom references in Fig. 2 and further information regarding RC determination is within the ESI† (page S9). US simulations used PLUMED V2.9.2's ASE interface.^23,24 Further details regarding simulation parameters can be found in the ESI† (page S11).

Table 1 Reaction coordinates for umbrella sampling

Reaction	RC definition	Sampled range (Å)
1A	d(15–6) + d(9–17)	6.50–2.40
1B	d(1–14) + d(4–6)	6.13–2.70
2	d((1–2) + (3–6)) + d(7–8)	7.12–3.71

---

Fig. 2 Atoms involved in transition state for each of reaction 1A, 1B and 2 are shown here in red.

The completion of US simulations along the RCs described in Table 1 resulted in good overlap of neighbouring windows resulting from WHAM analysis as shown in Fig. 3. Note the independence of each reaction's RC as shown in Table 1. As shown, sampling is performed well by ANI along the RC with a drop in sampling surrounding the high energy structures aligned with bond making and bond breaking along the RCs. This corresponds with the short lived low frequency nature of transition states, but highlights ANI-2x's ability to sample these high-energy infrequent conformers, albeit at reduced frequency compared to lower energy structures around the minima.

Fig. 3 Overlapping windows for each of (A) reaction 1A, (B) reaction 1B and (C) reaction 2, resulting from ANI-2x driven umbrella sampling. Plots show the conformer count within each window.

The resulting free energy surface (FES) from the biased US simulations for each reaction are shown in Fig. 4. The initial drop in energy at the start of each simulation is attributed to relaxation from slightly strained initial conformations.

Fig. 4 Free energy surfaces for (A) reaction 1A, (B) reaction 1B and (C) reaction 2 from ANI-2x driven umbrella sampling simulations.

Examination of the FES plots shown in Fig. 4, show clear differences from earlier static DFT results. Reaction 1A and 1B show increased barrier heights for the corresponding ANI-2x barriers. Reaction 1A performs well against the ANI-2x single point of the DFT transition state (62.1 versus 62.8 kcal mol⁻¹), however reaction 1B severely overestimates the transfer of the proton in the transition state (89.9 versus 33.3 kcal mol⁻¹). Reaction 2 shows an increased barrier compared to the ANI-2x single point of the transition state (52.9 versus 17.2 kcal mol⁻¹). Further, the FES show asynchronous transition states for each reaction where the DFT study above and resultant IRC calculations upon transition state structures have shown synchronous transition states.

The dynamic nature of US enables the discovery of new reaction pathways and transition states without requiring detailed prior knowledge of their exact conformations. While basic chemical intuition is needed to define the reaction coordinate, this approach allows us to systematically explore and identify potential transition states by sampling the relevant conformational space in an extensive manner. These initial transition state structures can then be further refined using high-level quantum mechanical calculations, offering a strategy for discovering novel reaction mechanisms without relying on precise initial guesses of transition state geometries or bond scanning calculations.

The energies predicted by US with ANI-2x differ from those of a static DFT approach, and are explored further within the ESI† (page S12). By plotting the bond distances along the reaction coordinate of each frame for the combined reaction 1A trajectory against the corresponding ANI-2x single-point energies (Fig. 5), we demonstrate thorough sampling of conformational space around the minima, followed by good sampling of high energy space surrounding the transition state.

Fig. 5 Energy of each frame in the US trajectory for reaction 1A, plotted against the reaction coordinate (Table 1).

Similar and additional analyses for reaction 1B and reaction 2 are shown in the ESI† (Fig. S7–S11). Finally, manual examination of the complete US trajectories for each reaction confirm that the reaction followed the expected RCs, successfully crossing the TS to reach the products. Furthermore, ANI-2x enables enhanced sampling that would be computationally prohibitive with DFT. Our US simulations explored 500 000 conformations across 50 windows per reaction in under 5 hours using 10 CPU cores, a scale made possible by ANI-2x's ∼27 000× speed advantage over DFT (0.02 vs. 552 seconds per single-point calculation).²

In conclusion, we have demonstrated the capability of ANI-2x for enhanced sampling in two distinct reactions, exploring three molecular transformations with varying conformational flexibility. While ANI-2x falls short of chemical accuracy versus DFT for high energy structures in this case, it excels in rapidly navigating the free energy surface and its ability to identify high-energy pathways provides valuable, systematically guided guesses for transition state calculations at higher levels of theory. Furthermore, ANI-2x effectively handles complex reaction mechanisms involving both synchronous and asynchronous bond rearrangements. Its rapid sampling of the FES and mechanistic insights highlight its utility for exploring conformationally flexible systems. Given that large biochemical systems, such as proteins and enzymes, are often intractable with pure DFT, ANI-2x's scalability and efficiency make it well suited for such applications.

Results were obtained using the ARCHIE-WeSt high performance computer (https://www.archie-west.ac.uk) based at the University of Strathclyde.

Ross Urquhart: methodology, validation, software, investigation, data curation, writing. Alexander van Teijlingen: methodology, validation, software, investigation, manuscript editing. Tell Tuttle: conceptualization, methodology, supervision and manuscript editing.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

Relevant data underpinning this work are available from the University of Strathclyde KnowledgeBase at: https://doi.org/10.15129/5c703b29-b0b84086-b261-af141d07e178. All code relevant to this publication can be accessed openly at https://github.com/Tuttlelab/ANI2xplorer.

Notes and references

1. J. S. Smith, O. Isayev and A. E. Roitberg, Chem. Sci., 2017, 8, 3192–3203.
2. C. Devereux, J. S. Smith, K. K. Huddleston, K. Barros, R. Zubatyuk, O. Isayev and A. E. Roitberg, J. Chem. Theory Comput., 2020, 16, 4192–4202.
3. I. Batatia, D. P. Kovacs, G. Simm, C. Ortner and G. Csanyi, Adv. Neural Inf. Process. Syst., 2022, 11423–11436.
4. R. Zubatyuk, J. S. Smith, J. Leszczynski and O. Isayev, Sci. Adv., 2019, 5, eaav6490.
5. R. J. Urquhart, A. Van Teijlingen and T. Tuttle, Phys. Chem. Chem. Phys., 2024, 26, 23934–23943.
6. J. Behler, J. Chem. Phys., 2016, 145, 170901.
7. S. Käser, L. I. Vazquez-Salazar, M. Meuwly and K. Töpfer, Digital Discovery, 2023, 2, 28–58.
8. H. Zhang, V. Juraskova and F. Duarte, Nat. Commun., 2024, 15, 6114.
9. L. L. Schaaf, E. Fako, S. De, A. Schäfer and G. Csányi, npj Comput. Mater., 2023, 9, 1–10.
10. S. Jiang, Y. Liu, C. Wang and T. Huang, Int. J. Quantum Chem., 2023, 123, e27087.
11. S. Zhang, R. Zubatyuk, Y. Yang, A. Roitberg and O. Isayev, J. Chem. Theory Comput., 2025, 21(9), 4365–4374.
12. J. Kästner, WIRES Comput. Mol. Sci., 2011, 1, 932–942.
13. J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2008, 128, 084106.
14. V. A. Rassolov, J. A. Pople, M. A. Ratner and T. L. Windus, J. Chem. Phys., 1998, 109, 1223–1229.
15. F. Neese, WIRES Comput. Mol. Sci., 2012, 2, 73–78.
16. F. Neese, WIRES Comput. Mol. Sci., 2022, 12, e1606.
17. K. Huddleston, R. Zubatyuk, J. Smith, A. Roitberg, O. Isayev, I. Pickering, C. Devereux and K. Barros, ANI-2x Release, 2023 DOI:10.5281/zenodo.10108942.
18. X. Gao, F. Ramezanghorbani, O. Isayev, J. S. Smith and A. E. Roitberg, J. Chem. Inf. Model., 2020, 60, 3408–3415.
19. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Duak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys.: Condens. Matter, 2017, 29, 273002.
20. E. C. B. Johnson and S. B. H. Kent, J. Am. Chem. Soc., 2006, 128, 6640–6646.
21. P. E. Dawson, T. W. Muir, I. Clark-Lewis and S. B. H. Kent, Science, 1994, 266, 776–779.
22. M. A. Hagras, M. A. Bellucci, G. Gobbo, R. A. Marek and B. L. Trout, J. Phys. Chem. B, 2020, 124, 9840–9851.
23. T. P. Consortium, Nat. Methods, 2019, 16, 670–673.
24. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comput. Phys. Commun., 2014, 185, 604–613.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5cc02090e

---