Enhanced Sampling for Free Energy Profiles with Post-Transition-State Bifurcations

Juno Nam YounJoon Jung [ [email protected]

Abstract

We present a method to explore the free energy landscapes of chemical reactions with post-transition-state bifurcations using an enhanced sampling method based on well-tempered metadynamics. Obviating the need for computationally expensive DFT-level ab initio molecular dynamics simulations, we obtain accurate energetics by utilizing a free energy perturbation scheme and deep learning estimator for the single-point energies of substrate configurations. Using a pair of easily interpretable collective variables, we present a quantitative free energy surface that is compatible with harmonic transition state theory calculations and in which the bifurcations are clearly visible. We demonstrate our approach with the example of the SpnF-catalyzed Diels–Alder reaction, a cycloaddition reaction in which post-transition-state bifurcation leads to the [4+2] as well as the [6+4] cycloadduct. We obtain the free energy landscapes for different stereochemical reaction pathways and characterize the mechanistic continuum between relevant reaction channels without explicitly searching for the pertinent transition state structures.

\hypersetup

colorlinks=true,allcolors=blue Seoul National University] Department of Chemistry, Seoul National University, Seoul 08826, South Korea

{tocentry} [Uncaptioned image]

Understanding the selectivity of a reaction is one of the most important challenges in computational chemistry. Reaction selectivity under kinetically controlled conditions can be predicted by estimating the relative heights of barriers or the free energies of transition states (TSs). Although this type of prediction is frequently utilized, its use implicitly assumes that competing reactive channels are well separated on the potential energy surface (PES), which is not always the case. Post-transition-state bifurcations (PTSBs) in chemical reactions occur when reactive trajectories passing through a single ambimodal TS bifurcate into two (or more) different products.^{1, 2, 3} Such phenomena have been observed for various organic, bioorganic, and organometallic systems.^{4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19}

The occurrence of bifurcations on PESs means that intrinsic reaction coordinate (IRC) analysis,^{20, 21} which is commonly used to study reactions, is not sufficient to characterize reactive systems featuring PTSBs. Accordingly, molecular dynamics (MD) simulations—in particular, quasiclassical direct dynamics and quantum mechanics/molecular mechanics (QM/MM) simulations—have been utilized to understand reaction mechanisms and selectivity associated with PTSBs.^{22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32} Although quasiclassical direct dynamics simulations account for product distributions and offer real-time trajectories from which mechanistic properties (e.g., synchronicity) can be probed, the results of such analyses are dependent on the starting TS structures. It has been reported that PESs in the vicinity of ambimodal transition states may possess a large number of saddle points because of the complexity of the chemical interactions and conformational flexibility of the reactants.³³ In this respect, consideration of the transition state “ensemble” is indispensable for obtaining a comprehensive understanding of complicated chemical reactions with PTSBs.

In this study, we aimed to explore free energy profiles in systems characterized by PTSBs by implicitly considering the relevant TS structures using an enhanced sampling method, well-tempered metadynamics (WTMetaD),³⁴ with ab initio molecular dynamics (AIMD).³⁵ Enhanced sampling methods based on metadynamics³⁶ have been increasingly utilized in recent years to explore the free energy landscapes of chemical reactions.^{37, 38, 39, 40, 41, 42} These methods allow efficient exploration of substrate configurations over the reaction barrier on the timescale of AIMD by accumulating a bias potential along predefined collective variables (CVs). Because of the high computational cost of AIMD, chemical reactions are often simulated at a semiempirical level, and hence the free energy surfaces (FESs) obtained from these dynamics simulations have intrinsically limited accuracy. Recently, the free energy perturbation (FEP) method has been applied to such simulations^{43, 44} to obtain more accurate FESs using single-point calculations for sampled configurations at a higher level of theory. Despite these improvements, as the number of atoms in the reactant increases, the increase in the computational cost of performing single-point calculations for a sufficiently large number of configurations to obtain a FES can be prohibitive. Motivated by the successful application of geometric deep learning methods for the prediction of various molecular properties,⁴⁵ we constructed a FES by training a graph neural network on the DFT calculation results for a small number of configurations sampled by WTMetaD simulation and predicting the energy for a large number of configurations.

To illustrate the application of our method to explore free energy landscapes with PTSBs, we characterized the ensemble of ambimodal transition states and relevant free energy landscapes of the SpnF-catalyzed Diels–Alder (DA) reaction (Figure 1a). The reaction corresponds to a [4+2] cycloaddition step in the biosynthesis of Spinosyn A and is characterized by PTSB because of the existence of ambimodal [6+4]/[4+2] bis-pericyclic (BPC) transition states.^{46, 47, 12} While several previous studies have focused on the effect of the enzyme and solvent on the reaction pathway and product selectivity (i.e., whether the reaction proceeds via BPC-type or DA-type TSs),^{22, 23, 24} Medvedev et al.³³ found that more than 300 BPC/DA TS structures can be located for the reaction, and the distinction between the BPC and DA mechanisms is blurred. Using our enhanced sampling method, we aimed to describe the “mechanistic continuum” by implicitly taking into account the relevant TS structures. By utilizing a pair of easily interpretable CVs, we obtained FESs in which the bifurcations were clearly visualized and consistent with previously located TS structures.

Refer to caption — Figure 1: SpnF-catalyzed Diels–Alder reaction. (a) Possible cycloaddition routes for the reaction. (b) Reactant conformations and [6+4] products with different stereochemistry. (c) Relevant geometric parameters (distances and torsional angles) for the substrate.

The effectiveness of a WTMetaD simulation depends substantially on the choice of CVs: they must act as low-dimensional descriptors that can discriminate the relevant metastable states and transition states.^{48, 49} Using the coordination numbers⁵⁰ given by

c_{i}=\frac{1-(d_{i}/\sigma)^{6}}{1-(d_{i}/\sigma)^{8}}\qquad(i\in\{1,2,3\})

(1)

where $d_{i}$ is the distance between two atoms as defined in Figure 1c and $\sigma$ is a reference distance (1.7 Å), we employed a set of two CVs, $\mathbf{s}=(s_{1},s_{2})$ , where $s_{1}=c_{1}+c_{2}+c_{3}$ follows the overall reaction progress and $s_{2}=c_{2}-c_{3}$ discriminates between the two products.

Although 3-A is considered to be the dominant [6+4] product¹² among the possible products with different stereochemistry (Figure 1b), an extensive TS search by Medvedev et al.³³ suggested the existence of several low-energy TS structures that are stereochemically relevant to other products (3-B–D). Accordingly, in our metadynamics simulations, we considered all the possible [6+4] products. However, since the transitions between metastable reactant states (1-A–D) leading to such products are hindered by the barriers between s-cis and s-trans isomers, using only the aforementioned CVs would require excessively long simulations.⁵¹ To overcome this problem, we introduced a set of conformational restraints $W_{\text{X}}(\mathbf{R})$ ( $\text{X}\in\{\mathrm{A,B,C,D}\}$ ) to constrain the substrate from accessing states other than 1-X, 2, and 3-X, and independent WTMetaD simulations with each of the restraints.

From the WTMetaD simulations, the FES can be reconstructed up to a constant $C_{\mathrm{X}}$ by reweighting⁵²:

	$\displaystyle F_{\mathrm{X}}(\mathbf{s})$	$\displaystyle=-\frac{1}{\beta}\log\left<\delta\left[\mathbf{s}-\mathbf{s}(\mathbf{R})\right]\right>+C_{\mathrm{X}}$
		$\displaystyle=-\frac{1}{\beta}\log\frac{\left<\delta\left[\mathbf{s}-\mathbf{s}(\mathbf{R})\right]e^{+\beta V(\mathbf{s}(\mathbf{R}))}\right>_{\mathrm{X}}}{\left<e^{+\beta V(\mathbf{s}(\mathbf{R}))}\right>_{\mathrm{X}}}+C_{\mathrm{X}}$		(2)

where $\left<\cdot\right>_{\mathrm{X}}$ denotes an average over the biased ensemble constructed from the WTMetaD run with restraint $W_{\mathrm{X}}(\mathbf{R})$ . To level off the FES from each WTMetaD run, we set the entire conformations of 1 as the reference state (i.e., $F_{\mathbf{1}}=0$ ). Then, we introduced another WTMetaD run (run T) with a set of three torsional angle CVs, $\Phi=(\phi_{1},\phi_{2},\phi_{3})$ (as defined in Figure 1c), to effectively sample the conformations of 1 over the s-cis/s-trans barriers. The free energy of state 1-X ( $\mathrm{X}\in\{\mathrm{A,B,C,D}\}$ ) can be calculated using the Zwanzig equation⁵³ as

e^{-\beta F_{\text{{1-X}}}}=\left<e^{-\beta W_{\mathrm{X}}(\mathbf{R})}\right>_{\mathbf{1}}=\frac{\left<e^{-\beta W_{\mathrm{X}}(\mathbf{R})}e^{+\beta V(\Phi(\mathbf{R}))}\right>_{\mathrm{T}}}{\left<e^{+\beta V(\Phi(\mathbf{R}))}\right>_{\mathrm{T}}}

(3)

where $\left<\cdot\right>_{\mathbf{1}}$ denotes an ensemble average over state 1 and $\left<\cdot\right>_{\mathrm{T}}$ a biased ensemble average obtained from WTMetaD run T. Finally, $C_{\mathrm{X}}$ in eq 2 can be determined by equating the free energy in eq 3 with the integral of the FES over $\Omega_{\text{{1-X}}}$ , the CV space corresponding to 1-X:

e^{-\beta F_{\text{{1-X}}}}=\int_{\Omega_{\text{{1-X}}}}e^{-\beta F_{\mathrm{X}}(\mathbf{s})}\,\mathrm{d}\mathbf{s}

(4)

Initially, 1.8-ns WTMetaD simulations were run using a semiempirical tight-binding Hamiltonian, GFN-xTB.⁵⁴ The computational cost of the semiempirical Hamiltonian was low enough to run sufficiently long simulations to sample the configuration space; however, the method does not provide a quantitatively accurate description of the chemical process of interest. Therefore, we employed the FEP scheme^{43, 44} to refine the FES at the DFT level of theory [M06-2X/6-31G(d)^{55, 56, 57, 58}]. The ensemble average of an observable $O(\mathbf{R})$ at a high-level Hamiltonian $U_{H}(\mathbf{R})$ can be calculated from ensemble averages using a low-level Hamiltonian $U_{L}(\mathbf{R})$ by incorporating the perturbative weight:

\left<O(\mathbf{R})\right>_{H}=\frac{\left<O(\mathbf{R})e^{+\beta[U_{L}(\mathbf{R})-U_{H}(\mathbf{R})]}\right>_{L}}{\left<e^{+\beta[U_{L}(\mathbf{R})-U_{H}(\mathbf{R})]}\right>_{L}}

(5)

where $\left<\cdot\right>_{H}$ and $\left<\cdot\right>_{L}$ are ensemble averages obtained using high-level and low-level Hamiltonians, respectively. Accordingly, the biased ensemble averages on the right-hand sides of eqs 2 and 3 are refined using eq 5.

Previous methods based on FEP^{44, 42} sampled a small number of configurations to calculate the DFT energies and perturbative weights and then applied post-processing such as moving average calculations to smooth the obtained high-level FES. As an alternative approach, in this study we utilized a graph neural network method inspired by molecular mechanics⁵⁹ trained on the DFT calculation results for a small subset of the configurations sampled from WTMetaD runs to estimate the DFT energies of a large number of configurations. The network was trained to simultaneously predict the differences between DFT and GFN-xTB energies, solvation energies in implicit water (CPCM^{60, 61}), and CM5 partial atomic charges⁶². As well as offering additional information for the system, such as the effect of implicit solvent, the multi-task setting⁶³ also provides more accurate DFT energy predictions (see Table S2). As shown in Figure 2a, the GFN-xTB and DFT energies are correlated but quantitatively different, and the values are also dependent on the reaction progress (CV $s_{1}$ ). The constructed deep learning model can accurately predict the DFT energy with a mean absolute error (MAE) of 0.486 kcal/mol on the test set (Figure 2b). Furthermore, the prediction error is less than 2 kcal/mol for more than 99% of the test set configurations, and the magnitude of the error does not depend on the reaction progress (see Figure S4). Hence, the deep learning model can be used as a reliable surrogate for the high-cost DFT single-point energy calculations, and when applied to FEP calculations, the resulting FES reflects the relatively high TS and product energies obtained from the DFT-level calculations (Figure 2c). Details of the restraints, WTMetaD simulations, and deep learning models are reported in the \hyperref[sec:methods]Computational Methods section.

Two-dimensional free energy profiles, $F(s_{1},s_{2})$ , obtained from each WTMetaD run (A, B, C, and D) with FEP refinement using eqs 2, 3, and 5 are shown in Figure 3. The free energy surfaces provide a quantitative description of the asymmetric post-TS bifurcations, and cross-sections along the reaction progress (CV $s_{1}$ ) illustrate the valley–ridge inflection point on the FES. The free energy landscapes near the transition region, especially the skewness relative to $s_{2}=0$ and the activation barriers, depend on the stereochemical pathways (A–D) determined by the respective stereochemical restraints. To better understand the product selectivity, we calculated the minimum free energy paths (MFEPs) using the nudged elastic band (NEB) algorithm.⁶⁴ MFEPs between the free energy minima of the reactants and products show that the pre-TS pathway for A aligns well with the $s_{2}=0$ line and the transition region has more “BPC-like” character, whereas the pathways for B–D are tilted toward the [4+2] product (2, $s_{2}>0$ ) and the transition region is “DA-like.” Although the MFEPs for B–D are similar, the curvatures of the FESs near the TS region are different, and a two-dimensional description is required to account for the product selectivity. We also obtained the TS structures at the current level of theory, by re-optimizing the TS structures reported by Medvedev et al.,³³ and projected them onto the CV space of the relevant stereochemistry (A–D). From Figure 3 and the TS energies given in Figure S3, it can be seen that free energy landscapes and associated MFEPs from the WTMetaD simulations (and FEP refinement) represent the energetic landscapes of the searched TS structures well. This result indicates that the free energy landscapes, reaction energetics, and preferred reaction pathways and mechanisms of the system considered were efficiently characterized using our methodology, obviating the need for an extensive manual search for TS structures.

In order to investigate the differences between the activation barriers of the different stereochemical pathways, one-dimensional profiles of the free energy versus the reaction progress were calculated (Figure 4). The results indicate that pathway A, with the lowest barrier, is the dominant stereochemical pathway, and the barrier increases in the order $\text{A}<\text{B}\ll\text{C}<\text{D}$ . As the barrier height depends on the choice of CVs, the gauge correction scheme described by Bal et al.⁶⁵ should be introduced for consistency. The free energy barrier $\Delta^{\ddagger}F$ is given as

\Delta^{\ddagger}F=F_{\text{TS}}^{G}+\frac{1}{\beta}\log\frac{\lambda}{h}\sqrt{\frac{2\pi m}{\beta}}

(6)

with the gauge-invariant geometric FES⁶⁶ $F^{G}(s_{1})$ defined as

F^{G}(s_{1})=F(s_{1})-\frac{1}{\beta}\log\left<\lambda\left|\nabla_{\mathbf{R}}s_{1}\right|\right>_{s_{1}}

(7)

where $\lambda$ is the length unit for the system coordinates $\mathbf{R}$ and $m$ is the reduced mass for the pair of atoms defining the CVs. The free energies of states 1-X (eq 3) and the free energy barriers $\Delta^{\ddagger}F_{\text{X}}$ (eq 6) for the respective stereochemical pathways ( $\text{X}\in\{\mathrm{A,B,C,D}\}$ ) are reported in Table 1.

Table 1: Free Energies of the Reactant States and Activation Free Energies for Different Stereochemical Pathways^a

Pathway (X)	$F_{\text{{1-X}}}$ (kcal/mol)	$\Delta^{\ddagger}F_{\text{X}}$ (kcal/mol)
A	$3.2\pm 0.6$	$28.3\pm 1.4$
B	$6.1\pm 1.1$	$29.6\pm 1.4$
C	$5.8\pm 0.8$	$32.1\pm 2.1$
D	$8.0\pm 1.8$	$33.2\pm 2.7$

^aThe free energies of the reactant states are reported with the standard deviation from five independent simulations. The uncertainty was propagated to the activation free energies since they were calculated from the FES shifted using $F_{\text{{1-X}}}$ (eq 4).

This result quantifies the aforementioned trend in the activation barriers of the different stereochemical pathways, with the lowest activation free energy, $\Delta^{\ddagger}F_{\text{A}}$ , being 28.3 kcal/mol. Only small changes in the reaction energetics were obtained when the solvation energy in implicit water was also applied at the FEP refining step (Table S1), with $\Delta^{\ddagger}F_{\text{A}}$ becoming 28.2 kcal/mol. Despite this activation barrier being similar to the value obtained using harmonic transition state theory and the same DFT functional (27.6 kcal/mol),¹² it is higher than the experimental activation barrier of 22.0 kcal/mol in water.⁴⁶ As the relative stability of configurations can be affected by hydrogen bonding,^{24, 23} QM/MM simulations are required for a more accurate description of the reaction.

Previous quasi-classical dynamics studies^{12, 24} have found that the [4+2]/[6+4] cycloadditions that are the focus of our study are energetically concerted but are characterized by stepwise dynamics. This has also been verified by means of kinetic isotope effect experiments.⁶⁷ To determine the step-wise bond formation dynamics, we constructed the FES with another set of CVs: $c_{1}$ and $\max(c_{2},c_{3})$ . Since these two CVs follow the formation of the first and second bonds during the cycloadditions, the resulting FES can be regarded as a More O’Ferrall–Jencks plot.⁶⁸ The obtained FES, the MFEP between the free energy minima, and the projected TS locations for stereochemical pathway A are shown in Figure 5. The reaction pathway, as represented by the MFEP, is curved toward the bottom right corner and indicates the two bonds are formed asynchronously. The free energy maximum along the reaction pathway coincides with the projected TS positions and corresponds to a substrate geometry in which the first bond is almost formed and the second is not. Similar results were obtained for the other stereochemical pathways (B–D), and these are reported in Figure S2.

We proposed a method to explore DFT-level free energy landscapes with post-TS bifurcation, using WTMetaD simulation and FEP refinement assisted by a deep learning energy estimator. Using our methodology, it is possible to characterize a bifurcating reaction pathway in a manner independent of a specific TS structure, as was demonstrated for the case of a SpnF-catalyzed Diels–Alder reaction. In addition, we showed that deep learning methods, especially graph neural networks, can be a reliable surrogate for costly single-point DFT calculations for FEP refining steps. In future studies, it would be advantageous to make use of specialized CVs^{69, 41, 70} and recently developed enhanced sampling methods^{71, 72} to study more complex cases (e.g., systems with post-TS trifurcations²⁹). Furthermore, to obtain a more accurate description of the reaction rates and product selectivity in these systems, future studies may benefit from the incorporation of explicit kinetics modeling.^{73, 74}

1 Computational Methods

Conformational restraints. While no restraint was applied to WTMetaD run T, two types of restraint were used for the other WTMetaD runs (A, B, C, D): one to prevent unwanted reactions, and the other to restrict the reactant to the stereochemical pathway of interest.

We defined the coordination number $c_{ij}\in[0,1]$ , angle $\theta_{ijkl}\in[0,\pi]$ , and torsional angle $\phi_{ijkl}\in(-\pi,\pi]$ as follows:

$\displaystyle c_{ij}$	$\displaystyle=\frac{1-(\|\mathbf{r}_{ij}\|/\sigma_{ij})^{6}}{1-(\|\mathbf{r}_{ij}\|/\sigma_{ij})^{8}}$	(8)
$\displaystyle\theta_{ijkl}$	$\displaystyle=\cos^{-1}\left(\frac{\mathbf{r}_{ji}\cdot\mathbf{r}_{kl}}{\left\|\mathbf{r}_{ji}\right\|\left\|\mathbf{r}_{kl}\right\|}\right)$	(9)
$\displaystyle\phi_{ijkl}$	$\displaystyle=\operatorname{atan2}\left(\left\|\mathbf{r}_{jk}\right\|\mathbf{r}_{ij}\cdot(\mathbf{r}_{jk}\times\mathbf{r}_{kl}),(\mathbf{r}_{ij}\times\mathbf{r}_{jk})\cdot(\mathbf{r}_{jk}\times\mathbf{r}_{kl})\right)$	(10)

where $\sigma_{ij}$ is 1.7 Å and 1.6 Å for the C–C and C–O atom pairs, respectively. First, to avoid the sampling of irrelevant reaction pathways, we restrained four atom pairs, $\mathcal{P}=\left\{(2,11),(3,14),(1,13),(11,19)\right\}$ , from being too close using the harmonic wall restraint

W_{\text{reac}}=k\sum_{(i,j)\in\mathcal{P}}\max(c_{ij}-0.4,\,0)^{2}

(11)

with $k=5\times 10^{3}$ kcal/mol. Further, to constrain the molecule to exploration of the relevant stereochemical pathway ( $\text{X}\in\{\mathrm{A,B,C,D}\}$ ), we also introduced the following harmonic wall restraint:

W_{\text{stereo,X}}=\kappa\sum_{i=1}^{10}\max(\alpha_{i,\text{X}}\cos\varphi_{i},\,0)^{2}

(12)

with $\kappa=50$ kcal/mol and the coefficients $\alpha_{i,\text{X}}$ and angles $\varphi_{i}$ defined as in Table 2; $\alpha_{i,\text{X}}$ values of $-1$ and $+1$ restrain the angle $\phi_{i}$ to values corresponding to the syn and anti configurations, respectively; $\varphi_{1}$ to $\varphi_{7}$ determine stereochemistry of the [4+2] product (2); and $\varphi_{8}$ to $\varphi_{10}$ control the stereochemical pathway leading to the [6+4] product (3).

$i$	$\varphi_{i}$	$\alpha_{i,\text{A}}$	$\alpha_{i,\text{B}}$	$\alpha_{i,\text{C}}$	$\alpha_{i,\text{D}}$
1	$\phi_{34,5,9,40}$	$-1$	$-1$	$-1$	$-1$
2	$\phi_{33,4,9,40}$	$+1$	$+1$	$+1$	$+1$
3	$\phi_{35,6,5,34}$	$+1$	$+1$	$+1$	$+1$
4	$\phi_{32,1,4,33}$	$-1$	$-1$	$-1$	$-1$
5	$\phi_{1,2,3,4}$	$-1$	$-1$	$-1$	$-1$
6	$\theta_{34,5,4,33}$	$+1$	$+1$	$+1$	$+1$
7	$\theta_{35,6,1,32}$	$-1$	$-1$	$-1$	$-1$
8	$\phi_{43,14,6,35}$	$-1$	$-1$	$+1$	$+1$
9	$\phi_{44,18,1,32}$	$-1$	$+1$	$-1$	$+1$
10	$\theta_{43,14,18,44}$	$-1$	$+1$	$+1$	$-1$

Then, the restraints for the WTMetaD runs, mentioned in the main text, are defined by $W_{\text{X}}(\mathbf{R})=W_{\text{reac}}+W_{\text{stereo,X}}$ . Although the stereochemical restraint may prevent the sampling of the full conformational space of the products, its effect on the obtained FES in this work is insignificant, and this is apparent from the consistency of the FESs in the product region in Figure 4.

Simulation settings. WTMetaD simulations were carried out using the CP2K 9.1 software package⁷⁵ patched with PLUMED 2.7.3.⁷⁶ The semiempirical GFN-xTB Hamiltonian⁵⁴ was used to obtain trajectories with an integration time step of 0.5 fs. A canonical ( $NVT$ ) ensemble at 298.15 K was sampled using the velocity rescaling thermostat⁷⁷ with a time constant of 50 fs. For each WTMetaD run (A, B, C, D, or T) type, the initial structure was obtained by energy minimization at the GFN-xTB level. The system was equilibrated under the application of the respective constraint $W_{\text{X}}(\mathbf{R})$ for 1.0 ns, and five configurations (positions and velocities) were sampled from the last 0.8 ns to initiate five statistically independent WTMetaD runs for each type. For WTMetaD runs A, B, C, and D, a bias factor ( $\gamma$ ) of 30 was used, and Gaussian biases with a height of 1.0 kcal/mol and width of 0.05 were added every 200 MD steps. For WTMetaD run T, a bias factor of 15 was used, and Gaussian biases with a height of 0.5 kcal/mol and a width of 0.2 rad were added every 200 MD steps. All the WTMetaD simulations were run for 1.8 ns, and the last 1.2 ns was used to obtain the free energy profiles.

DFT calculations. The single point energy calculations for the sampled configurations and the TS structure optimizations were carried out at the M06-2X/6-31G(d) level of theory^{55, 56, 57, 58} as implemented in Gaussian 16, revision C.01.⁷⁸ The keywords used to designate the level of theory and implicit solvation were m062x/6-31g(d) and scrf=(cpcm,solvent=water), respectively. We used the keyword pop=hirshfeld to obtain CM5 partial charges. The transition states were optimized using the keyword opt=(ts,calcfc,noeigen).

The transition state structures in ref 33 were obtained at the M06-2X/6-31+G(d) level of theory. We re-optimized these transition state structures at our level of theory [M06-2X/6-31G(d)] to verify that the transition regions in the free energy landscapes obtained in our work were consistent with the results of harmonic transition state theory. With the lowest energy TS as a reference, only TSs with energies of less than 15 kcal/mol were used. The projected TS locations are shown in the upper panels of Figure 3, and the relative energies are reported in Figure S3.

Deep learning model. The geometric deep learning model used in this work was MXMNet.⁵⁹ While other similar models (e.g., DimeNet⁺⁺⁷⁹) might have been considered suitable for our task, we did not benchmark the performance of other models as DFT-level energy estimation by deep learning was not the main source of error for the free energy surfaces obtained in this work.

The graph neural network used in this study simultaneously predicted (1) the differences between DFT and GFN-xTB energies ( $E_{\text{DFT}}-E_{\text{xTB}}$ ), (2) the solvation energies in implicit water ( $E_{\text{DFT/solv}}-E_{\text{DFT}}$ ), and (3) the CM5 partial atomic charges ( $\mathbf{q}$ ). Accordingly, the MXMNet architecture was modified so that each Local Layer Message Passing Module had three separate “Output” layers, one for each property to be predicted (Figure 3c in ref 59). The target values for the energies were obtained by subtracting the reference values obtained from the molecular geometry with the lowest $E_{\text{DFT}}$ . Energy value units of kcal/mol were used, along with partial charge units of $e$ .

The deep learning model was optimized by minimizing the weighted sum of the L1 losses for the targets, i.e.,

\mathcal{L}=L_{1}(E_{\text{DFT}}-E_{\text{xTB}})+w_{s}\cdot L_{1}(E_{\text{DFT/solv}}-E_{\text{DFT}})+w_{c}\cdot L_{1}(\mathbf{q})

(13)

where $L_{1}(\cdot)$ denotes the mean absolute error between the predicted and true values. Several values for the weights $w_{s}$ and $w_{c}$ were tested as in Table S2, and the weights for the final ensemble were determined as $w_{s}=1$ and $w_{c}=2500$ . The loss was optimized by the Adam optimizer⁸⁰ with an initial learning rate of $3\times 10^{-4}$ . The learning rate was warmed up linearly from zero to the initial learning rate for the first five epochs, then decayed exponentially by a factor of $0.5^{1/250}$ every epoch. We used a batch size of 64 and a hidden dimension size of 128; the number of layers was 6, the gradient clipping norm was 100, and the global and local layer distance cutoff values were 10 Å and 5 Å, respectively.

{acknowledgement}

The authors thank Jaehak Lee for their careful reading of the manuscript and constructive suggestions.

{suppinfo}

Free energy profiles calculated at the GFN-xTB level using different stereochemical restraints, with corresponding transition state energies; calculated free energies with implicit solvation; error distribution of the deep learning model outputs

References

Ess et al. 2008 Ess, D. H.; Wheeler, S. E.; Iafe, R. G.; Xu, L.; Çelebi-Ölçüm, N.; Houk, K. N. Bifurcations on Potential Energy Surfaces of Organic Reactions. Angew. Chem., Int. Ed. 2008, 47, 7592–7601
Rehbein and Carpenter 2011 Rehbein, J.; Carpenter, B. K. Do We Fully Understand What Controls Chemical Selectivity? Phys. Chem. Chem. Phys. 2011, 13, 20906–20922
Hare and Tantillo 2017 Hare, S. R.; Tantillo, D. J. Post-Transition State Bifurcations Gain Momentum – Current State of the Field. Pure Appl. Chem. 2017, 89, 679–698
Yamataka et al. 1999 Yamataka, H.; Aida, M.; Dupuis, M. One Transition State Leading to Two Product States: Ab Initio Molecular Dynamics Simulations of the Reaction of Formaldehyde Radical Anion and Methyl Chloride. Chem. Phys. Lett. 1999, 300, 583–587
Caramella et al. 2002 Caramella, P.; Quadrelli, P.; Toma, L. An Unexpected Bispericyclic Transition Structure Leading to 4+2 and 2+4 Cycloadducts in the Endo Dimerization of Cyclopentadiene. J. Am. Chem. Soc. 2002, 124, 1130–1131
Ussing et al. 2006 Ussing, B. R.; Hang, C.; Singleton, D. A. Dynamic Effects on the Periselectivity, Rate, Isotope Effects, and Mechanism of Cycloadditions of Ketenes with Cyclopentadiene. J. Am. Chem. Soc. 2006, 128, 7594–7607
Çelebi-Ölçüm et al. 2007 Çelebi-Ölçüm, N.; Ess, D. H.; Aviyente, V.; Houk, K. Lewis Acid Catalysis Alters the Shapes and Products of Bis-Pericyclic Diels–Alder Transition States. J. Am. Chem. Soc. 2007, 129, 4528–4529
Thomas et al. 2008 Thomas, J. B.; Waas, J. R.; Harmata, M.; Singleton, D. A. Control Elements in Dynamically Determined Selectivity on a Bifurcating Surface. J. Am. Chem. Soc. 2008, 130, 14544–14555
Hong and Tantillo 2009 Hong, Y. J.; Tantillo, D. J. A Potential Energy Surface Bifurcation in Terpene Biosynthesis. Nat. Chem. 2009, 1, 384–389
Siebert et al. 2011 Siebert, M. R.; Zhang, J.; Addepalli, S. V.; Tantillo, D. J.; Hase, W. L. The Need for Enzymatic Steering in Abietic Acid Biosynthesis: Gas-Phase Chemical Dynamics Simulations of Carbocation Rearrangements on a Bifurcating Potential Energy Surface. J. Am. Chem. Soc. 2011, 133, 8335–8343
Zhang et al. 2015 Zhang, L.; Wang, Y.; Yao, Z.-J.; Wang, S.; Yu, Z.-X. Kinetic or Dynamic Control on a Bifurcating Potential Energy Surface? An Experimental and DFT Study of Gold-Catalyzed Ring Expansion and Spirocyclization of 2-Propargyl- $\beta$ -Tetrahydrocarbolines. J. Am. Chem. Soc. 2015, 137, 13290–13300
Patel et al. 2016 Patel, A.; Chen, Z.; Yang, Z.; Gutiérrez, O.; Liu, H.-w.; Houk, K.; Singleton, D. A. Dynamically Complex [6+4] and [4+2] Cycloadditions in the Biosynthesis of Spinosyn A. J. Am. Chem. Soc. 2016, 138, 3631–3634
Hare and Tantillo 2017 Hare, S. R.; Tantillo, D. J. Cryptic Post-Transition State Bifurcations that Reduce the Efficiency of Lactone-Forming Rh-Carbenoid C–H Insertions. Chem. Sci. 2017, 8, 1442–1449
Hare et al. 2017 Hare, S. R.; Pemberton, R. P.; Tantillo, D. J. Navigating Past a Fork in the Road: Carbocation- $\pi$ Interactions Can Manipulate Dynamic Behavior of Reactions Facing Post-Transition-State Bifurcations. J. Am. Chem. Soc. 2017, 139, 7485–7493
Yu et al. 2017 Yu, P.; Chen, T. Q.; Yang, Z.; He, C. Q.; Patel, A.; Lam, Y.-h.; Liu, C.-Y.; Houk, K. Mechanisms and Origins of Periselectivity of the Ambimodal [6+4] Cycloadditions of Tropone to Dimethylfulvene. J. Am. Chem. Soc. 2017, 139, 8251–8258
Hare et al. 2018 Hare, S. R.; Li, A.; Tantillo, D. J. Post-Transition State Bifurcations Induce Dynamical Detours in Pummerer-Like Reactions. Chem. Sci. 2018, 9, 8937–8945
Xue et al. 2019 Xue, X.-S.; Jamieson, C. S.; Garcia-Borràs, M.; Dong, X.; Yang, Z.; Houk, K. Ambimodal Trispericyclic Transition State and Dynamic Control of Periselectivity. J. Am. Chem. Soc. 2019, 141, 1217–1221
Liu et al. 2020 Liu, F.; Chen, Y.; Houk, K. Huisgen’s 1,3-Dipolar Cycloadditions to Fulvenes Proceed via Ambimodal [6+4]/[4+2] Transition States. Angew. Chem., Int. Ed. 2020, 59, 12412–12416
Zhang et al. 2021 Zhang, H.; Novak, A. J.; Jamieson, C. S.; Xue, X.-S.; Chen, S.; Trauner, D.; Houk, K. Computational Exploration of the Mechanism of Critical Steps in the Biomimetic Synthesis of Preuisolactone A, and Discovery of New Ambimodal (5+2)/(4+2) Cycloadditions. J. Am. Chem. Soc. 2021, 143, 6601–6608
Fukui 1970 Fukui, K. Formulation of the Reaction Coordinate. J. Phys. Chem. 1970, 74, 4161–4163
Fukui 1981 Fukui, K. The Path of Chemical Reactions — The IRC Approach. Acc. Chem. Res. 1981, 14, 363–368
Zheng and Thiel 2017 Zheng, Y.; Thiel, W. Computational Insights into an Enzyme-Catalyzed [4+2] Cycloaddition. J. Org. Chem. 2017, 82, 13563–13571
Chen et al. 2018 Chen, N.; Zhang, F.; Wu, R.; Hess Jr, B. A. Biosynthesis of Spinosyn A: A [4+2] or [6+4] Cycloaddition? ACS Catal. 2018, 8, 2353–2358
Yang et al. 2018 Yang, Z.; Yang, S.; Yu, P.; Li, Y.; Doubleday, C.; Park, J.; Patel, A.; Jeon, B.-s.; Russell, W. K.; Liu, H.-w. et al. Influence of Water and Enzyme SpnF on the Dynamics and Energetics of the Ambimodal [6+4]/[4+2] Cycloaddition. Proc. Natl. Acad. Sci. U. S. A. 2018, 115, E848–E855
Yang et al. 2018 Yang, Z.; Zou, L.; Yu, Y.; Liu, F.; Dong, X.; Houk, K. Molecular Dynamics of the Two-Stage Mechanism of Cyclopentadiene Dimerization: Concerted or Stepwise? Chem. Phys. 2018, 514, 120–125
Yang et al. 2018 Yang, Z.; Dong, X.; Yu, Y.; Yu, P.; Li, Y.; Jamieson, C.; Houk, K. Relationships Between Product Ratios in Ambimodal Pericyclic Reactions and Bond Lengths in Transition Structures. J. Am. Chem. Soc. 2018, 140, 3061–3067
Yang et al. 2019 Yang, Z.; Jamieson, C. S.; Xue, X.-S.; Garcia-Borràs, M.; Benton, T.; Dong, X.; Liu, F.; Houk, K. Mechanisms and Dynamics of Reactions Involving Entropic Intermediates. Trends Chem. 2019, 1, 22–34
Zou and Houk 2021 Zou, Y.; Houk, K. Mechanisms and Dynamics of Synthetic and Biosynthetic Formation of Delitschiapyrones: Solvent Control of Ambimodal Periselectivity. J. Am. Chem. Soc. 2021, 143, 11734–11740
Jamieson et al. 2021 Jamieson, C. S.; Sengupta, A.; Houk, K. Cycloadditions of Cyclopentadiene and Cycloheptatriene with Tropones: All Endo-[6+4] Cycloadditions are Ambimodal. J. Am. Chem. Soc. 2021, 143, 3918–3926
Wang et al. 2021 Wang, X.; Zhang, C.; Jiang, Y.; Wang, W.; Zhou, Y.; Chen, Y.; Zhang, B.; Tan, R. X.; Ge, H. M.; Yang, Z. J. et al. Influence of Water and Enzyme on the Post-Transition State Bifurcation of NgnD-Catalyzed Ambimodal [6+4]/[4+2] Cycloaddition. J. Am. Chem. Soc. 2021, 143, 21003–21009
Qin et al. 2021 Qin, Z.-X.; Tremblay, M.; Hong, X.; Yang, Z. J. Entropic Path Sampling: Computational Protocol to Evaluate Entropic Profile along a Reaction Path. J. Phys. Chem. Lett. 2021, 12, 10713–10719
Ito et al. 2022 Ito, T.; Maeda, S.; Harabuchi, Y. Kinetic Analysis of a Reaction Path Network Including Ambimodal Transition States: A Case Study of an Intramolecular Diels–Alder Reaction. J. Chem. Theory Comput. 2022, 18, 1663–1671
Medvedev et al. 2017 Medvedev, M. G.; Zeifman, A. A.; Novikov, F. N.; Bushmarinov, I. S.; Stroganov, O. V.; Titov, I. Y.; Chilov, G. G.; Svitanko, I. V. Quantifying Possible Routes for SpnF-Catalyzed Formal Diels–Alder Cycloaddition. J. Am. Chem. Soc. 2017, 139, 3942–3945
Barducci et al. 2008 Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 020603
Marx and Hutter 2009 Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods; Cambridge University Press: Cambridge, U.K., 2009
Laio and Parrinello 2002 Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562–12566
Piccini et al. 2018 Piccini, G.; Mendels, D.; Parrinello, M. Metadynamics with Discriminants: A Tool for Understanding Chemistry. J. Chem. Theory Comput. 2018, 14, 5040–5044
Rizzi et al. 2019 Rizzi, V.; Mendels, D.; Sicilia, E.; Parrinello, M. Blind Search for Complex Chemical Pathways Using Harmonic Linear Discriminant Analysis. J. Chem. Theory Comput. 2019, 15, 4507–4515
Schilling et al. 2020 Schilling, M.; Cunha, R. A.; Luber, S. Zooming in on the O–O Bond Formation—–An Ab Initio Molecular Dynamics Study Applying Enhanced Sampling Techniques. J. Chem. Theory Comput. 2020, 16, 2436–2449
Fu et al. 2021 Fu, Y.; Bernasconi, L.; Liu, P. Ab Initio Molecular Dynamics Simulations of the S_N1/S_N2 Mechanistic Continuum in Glycosylation Reactions. J. Am. Chem. Soc. 2021, 143, 1577–1589
Trizio and Parrinello 2021 Trizio, E.; Parrinello, M. From Enhanced Sampling to Reaction Profiles. J. Phys. Chem. Lett. 2021, 12, 8621–8626
Raucci et al. 2022 Raucci, U.; Rizzi, V.; Parrinello, M. Discover, Sample, and Refine: Exploring Chemistry with Enhanced Sampling Techniques. J. Phys. Chem. Lett. 2022, 13, 1424–1430
Li et al. 2018 Li, P.; Jia, X.; Pan, X.; Shao, Y.; Mei, Y. Accelerated Computation of Free Energy Profile at Ab Initio Quantum Mechanical/Molecular Mechanics Accuracy via a Semi-Empirical Reference Potential. I. Weighted Thermodynamics Perturbation. J. Chem. Theory Comput. 2018, 14, 5583–5596
Piccini and Parrinello 2019 Piccini, G.; Parrinello, M. Accurate Quantum Chemical Free Energies at Affordable Cost. J. Phys. Chem. Lett. 2019, 10, 3727–3731
Atz et al. 2021 Atz, K.; Grisoni, F.; Schneider, G. Geometric Deep Learning on Molecular Representations. Nat. Mach. Intell. 2021, 3, 1023–1032
Kim et al. 2011 Kim, H. J.; Ruszczycky, M. W.; Choi, S.-h.; Liu, Y.-n.; Liu, H.-w. An Enzyme-Catalyzed [4+2] Cycloaddition is a Key Step in the Biosynthesis of Spinosyn A. Nature 2011, 473, 109–112
Fage et al. 2015 Fage, C. D.; Isiorho, E. A.; Liu, Y.; Wagner, D. T.; Liu, H.-w.; Keatinge-Clay, A. T. The Structure of SpnF, a Standalone Enzyme that Catalyzes [4+2] Cycloaddition. Nat. Chem. Biol. 2015, 11, 256–258
Valsson et al. 2016 Valsson, O.; Tiwary, P.; Parrinello, M. Enhancing Important Fluctuations: Rare Events and Metadynamics from a Conceptual Viewpoint. Annu. Rev. Phys. Chem. 2016, 67, 159–184
Bussi and Laio 2020 Bussi, G.; Laio, A. Using Metadynamics to Explore Complex Free-Energy Landscapes. Nat. Rev. Phys. 2020, 2, 200–212
Iannuzzi et al. 2003 Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2003, 90, 238302
Pietrucci 2017 Pietrucci, F. Strategies for the exploration of free energy landscapes: Unity in diversity and challenges ahead. Rev. Phys. 2017, 2, 32–45
Branduardi et al. 2012 Branduardi, D.; Bussi, G.; Parrinello, M. Metadynamics with Adaptive Gaussians. J. Chem. Theory Comput. 2012, 8, 2247–2254
Zwanzig 1954 Zwanzig, R. W. High-temperature equation of state by a perturbation method. I. Nonpolar gases. J. Chem. Phys. 1954, 22, 1420–1426
Grimme et al. 2017 Grimme, S.; Bannwarth, C.; Shushkov, P. A Robust and Accurate Tight-Binding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All spd-Block Elements ( $Z$ = 1–86). J. Chem. Theory Comput. 2017, 13, 1989–2009
Zhao and Truhlar 2008 Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215–241
Ditchfield et al. 1971 Ditchfield, R.; Hehre, W. J.; Pople, J. A. Self-Consistent Molecular-Orbital Methods. IX. An Extended Gaussian-Type Basis for Molecular-Orbital Studies of Organic Molecules. J. Chem. Phys. 1971, 54, 724–728
Hehre et al. 1972 Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self-Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. J. Chem. Phys. 1972, 56, 2257–2261
Hariharan and Pople 1973 Hariharan, P. C.; Pople, J. A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theor. Chim. Acta 1973, 28, 213–222
Zhang et al. 2020 Zhang, S.; Liu, Y.; Xie, L. Molecular Mechanics-Driven Graph Neural Network with Multiplex Graph for Molecular Structures. arXiv preprint, arXiv:2011.07457, 2020
Barone and Cossi 1998 Barone, V.; Cossi, M. Quantum Calculation of Molecular Energies and Energy Gradients in Solution by a Conductor Solvent Model. J. Phys. Chem. A 1998, 102, 1995–2001
Cossi et al. 2003 Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. Energies, Structures, and Electronic Properties of Molecules in Solution With the C-PCM Solvation Model. J. Comput. Chem. 2003, 24, 669–681
Marenich et al. 2012 Marenich, A. V.; Jerome, S. V.; Cramer, C. J.; Truhlar, D. G. Charge Model 5: An Extension of Hirshfeld Population Analysis for the Accurate Description of Molecular Interactions in Gaseous and Condensed Phases. J. Chem. Theory Comput. 2012, 8, 527–541
Caruana 1997 Caruana, R. Multitask Learning. Mach. Learn. 1997, 28, 41–75
Henkelman and Jónsson 2000 Henkelman, G.; Jónsson, H. Improved Tangent Estimate in the Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points. J. Chem. Phys. 2000, 113, 9978–9985
Bal et al. 2020 Bal, K. M.; Fukuhara, S.; Shibuta, Y.; Neyts, E. C. Free Energy Barriers from Biased Molecular Dynamics Simulations. J. Chem. Phys. 2020, 153, 114118
Hartmann and Schütte 2007 Hartmann, C.; Schütte, C. Comment on Two Distinct Notions of Free Energy. Phys. D 2007, 228, 59–63
Jeon et al. 2017 Jeon, B.-s.; Ruszczycky, M. W.; Russell, W. K.; Lin, G.-M.; Kim, N.; Choi, S.-h.; Wang, S.-A.; Liu, Y.-n.; Patrick, J. W.; Russell, D. H. et al. Investigation of the Mechanism of the SpnF-Catalyzed [4+2]-Cycloaddition Reaction in the Biosynthesis of Spinosyn A. Proc. Natl. Acad. Sci. U. S. A. 2017, 114, 10408–10413
Anslyn and Dougherty 2006 Anslyn, E. V.; Dougherty, D. A. Modern Physical Organic Chemistry; University Science Books: Sausalito, CA, 2006; pp 407–412
Bonati et al. 2020 Bonati, L.; Rizzi, V.; Parrinello, M. Data-Driven Collective Variables for Enhanced Sampling. J. Phys. Chem. Lett. 2020, 11, 2998–3004
Bonati et al. 2021 Bonati, L.; Piccini, G.; Parrinello, M. Deep Learning the Slow Modes for Rare Events Sampling. Proc. Natl. Acad. Sci. U. S. A. 2021, 118, e2113533118
Invernizzi and Parrinello 2020 Invernizzi, M.; Parrinello, M. Rethinking Metadynamics: From Bias Potentials to Probability Distributions. J. Phys. Chem. Lett. 2020, 11, 2731–2736
Invernizzi and Parrinello 2022 Invernizzi, M.; Parrinello, M. Exploration vs Convergence Speed in Adaptive-bias Enhanced Sampling. arXiv preprint, arXiv:2201.09950, 2022
Tiwary and Parrinello 2013 Tiwary, P.; Parrinello, M. From Metadynamics to Dynamics. Phys. Rev. Lett. 2013, 111, 230602
McCarty et al. 2015 McCarty, J.; Valsson, O.; Tiwary, P.; Parrinello, M. Variationally Optimized Free-Energy Flooding for Rate Calculation. Phys. Rev. Lett. 2015, 115, 070601
Kühne et al. 2020 Kühne, T. D.; Iannuzzi, M.; Del Ben, M.; Rybkin, V. V.; Seewald, P.; Stein, F.; Laino, T.; Khaliullin, R. Z.; Schütt, O.; Schiffmann, F. et al. CP2K: An Electronic Structure and Molecular Dynamics Software Package - Quickstep: Efficient and Accurate Electronic Structure Calculations. J. Chem. Phys. 2020, 152, 194103
Tribello et al. 2014 Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New Feathers for an Old Bird. Comput. Phys. Commun. 2014, 185, 604–613
Bussi et al. 2007 Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101
Frisch et al. 2016 Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H. et al. Gaussian 16, revision C.01. Gaussian, Inc: Wallingford, CT, 2016
Klicpera et al. 2020 Klicpera, J.; Giri, S.; Margraf, J. T.; Günnemann, S. Fast and Uncertainty-Aware Directional Message Passing for Non-Equilibrium Molecules. arXiv preprint, arXiv:2011.14115, 2020
Kingma and Ba 2014 Kingma, D. P.; Ba, J. Adam: A Method for Stochastic Optimization. arXiv preprint, arXiv:1412.6980, 2014