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

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.


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 welltempered 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.

Graphical TOC Entry
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) simulationsin 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. 33 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), 34 with ab initio molecular dynamics (AIMD). 35 Enhanced sampling methods based on metadynamics 36 have been increasingly utilized in recent years to explore the free energy landscapes of chemical reactions. [37][38][39][40][41][42] These meth-ods 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 states. 12,46,47 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. 33 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. 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 50 given by where d i is the distance between two atoms as defined in Figure 1c and σ is a reference distance (1.7 Å), we employed a set of two CVs, 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. simulations. 51 To overcome this problem, we introduced a set of conformational restraints W X (R) (X ∈ {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 X by reweighting: 52 where · X denotes an average over the biased ensemble constructed from the WTMetaD run with restraint W X (R). To level off the FES from each WTMetaD run, we set the entire conformations of 1 as the reference state (i.e., F 1 = 0). Then, we introduced another WTMetaD run (run T) with a set of three torsional angle CVs, Φ = (φ 1 , φ 2 , φ 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 (X ∈ {A, B, C, D}) can be calculated using the Zwanzig equation 53 where · 1 denotes an ensemble average over state 1 and · T a biased ensemble average obtained from WTMetaD run T. Finally, C X in eq 2 can be determined by equating the free energy in eq 3 with the integral of the FES over Ω 1-X , the CV space corresponding to 1-X: Initially, 1.8-ns WTMetaD simulations were run using a semiempirical tight-binding Hamiltonian, GFN-xTB. 54 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(R) at a high-level Hamiltonian U H (R) can be calculated from ensemble averages using a low-level Hamiltonian U L (R) by incorporating the perturbative weight: where · H and · 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 42,44 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 59 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. 62 As well as offering additional information for the system, such as the effect of implicit solvent, the multi-task setting 63 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 with the gauge-invariant geometric FES 66 F G (s 1 ) defined as where λ is the length unit for the system coordinates 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 ∆ ‡ F X (eq 6) for the respective stereochemical pathways (X ∈ {A, B, C, D}) are reported in Table 1. a The 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 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, ∆ ‡ F 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) 67 To determine the step-wise bond formation dynamics, we constructed the FES with  Figure S2. The formation of the two bonds is represented by the coordination number CVs, c 1 and max(c 2 , c 3 ) (see Figure 1c and eq 1). The minimum free energy path (MFEP) between the free energy minima and the transition state structures are displayed in the same way as those in Figure 3.  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 41,69,70 and recently developed enhanced sampling methods 71,72 to study more complex cases (e.g., systems with post-TS trifurcations 29 ). 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 Computational Methods We defined the coordination number c ij ∈ [0, 1], angle θ ijkl ∈ [0, π], and torsional angle φ ijkl ∈ (−π, π] as follows: θ ijkl = cos −1 r ji · r kl |r ji | |r kl | (9) where σ 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, P = {(2, 11), (3,14), (1,13), (11,19)}, from being too close using the harmonic wall restraint with k = 5 × 10 3 kcal/mol. Further, to constrain the molecule to exploration of the relevant stereochemical pathway (X ∈ {A, B, C, D}), we also introduced the following harmonic wall restraint: with κ = 50 kcal/mol and the coefficients α i,X and angles ϕ i defined as in Table 2  Then, the restraints for the WTMetaD runs, mentioned in the main text, are defined by W X (R) = W reac + W 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. The deep learning model was optimized by minimizing the weighted sum of the L1 losses for the targets, i.e., where L 1 (·) 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 80 with an initial learning rate of 3×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.   The absolute prediction error is less than 2 kcal/mol for more than 99% of the test set configurations. S-3