Toward Understanding CB[7]-Based Supramolecular Diels-Alder Catalysis

Cucurbiturils (CBs) are robust and versatile macrocyclic compounds, often used as molecular hosts in complex supramolecular systems. In previous work, remarkable catalytic activity has been observed for asymmetric cycloadditions under very mild conditions. Herein, we investigate the nature of supramolecular catalysis using DFT calculations and QM/MM techniques. We discuss induced conformational changes, electrostatic shielding effects from the highly polar aqueous environment and cooperativity in hydrogen bonding of the substrates in explicit water using QM/MM simulation techniques. Our results show little specificity for the chosen molecules, suggesting an excellent opportunity to expand the scope for catalytic use of these supramolecular macrocyclic containers.


INTRODUCTION
Sustainable catalysis is evolving to limit the environmental impact of chemical industries. Organocatalytic and supramolecular systems are successfully developed offering mild reaction conditions and great selectivity. These are often aimed at mimicking enzymes that achieve the highest catalytic efficiencies as well as specificity. Cucurbit[n]urils (CB[n], n = 5-8), are great supramolecular hosts, creating a microscopic environment while maintaining good solubility in water (Barrow et al., 2015). Besides participating in host-guest chemistry (Zhang et al., 2012;Wu et al., 2017Wu et al., , 2019Olesińska et al., 2019;Huang et al., 2020), they were successfully utilized in catalytic applications (Wang et al., 2011;Taylor et al., 2013;Zheng et al., 2015;Palma et al., 2017;Ren et al., 2018;Rad et al., 2019), surface chemistry (Wagner et al., 2020), and sensing and detection (Huang et al., 2012;Kasera et al., 2012;Biedermann et al., 2015;de Nijs et al., 2017de Nijs et al., , 2019. Substrate binding within CB[n]-cavities is usually observed spectroscopically (through a variety of NMR techniques, surface Raman spectroscopy as well as UV/Vis and PL spectroscopies) (del Barrio et al., 2016;Sigle et al., 2016) and using isothermal titration calorimetry (ITC) experiments. However, thermodynamic parameters and binding constants attainable by these experiments do not alone explain the catalytic effect and more structural insights are needed to further explore these phenomena. Computational studies can account for many aspects of the catalysis (Harvey et al., 2019) enhancing our understanding of the origins of these effects. Cucurbiturils are subject to many computational works, often in conjecture with experiments, contributing to interpretation of spectroscopy (Chio et al., 2019), electronic properties (Sin et al., 2020) and in particular binding free energies and structures (Yin et al., 2016;Senthilnathan et al., 2018;Olesińska et al., 2019;Wu et al., 2019). As hosts, CBs were featured in many blind prediction challenges (Muddana et al., 2012(Muddana et al., , 2014Lee et al., 2017;Rizzi et al., 2018).
Herein we focus on the CB-modulated catalysis of intramolecular Diels-Alder cycloaddition reactions (Palma et al., 2017) (Figure 1). N-allyl-2-furfurylamines are known model substrates for Diels-Alder reactions (Andrés et al., 1997;Gschwend et al., 2002). Based on the cavity size and our previously reported binding studies, CB[7] was found to be the best host homolog to offer enhanced water solubility allowing for higher concentrations, conformational restriction and tight binding. In physiological conditions, low conversion is achieved, while in the presence of CB[7], the cycloaddition rates are significantly sped up, and the reaction goes to completion within minutes. Importantly, at 10 mol% catalyst loading, the reaction still occurred at an increased rate, indicating a catalytic turnover in the order of organocatalysts. This work was the first example of carrying out this reaction in aqueous environment without incorporation of protecting groups and was compared to the activity of natural Diels-Alderase (Ose et al., 2003;Byrne et al., 2016).
We have previously found that binding in CB favors the closed conformation of the reactants, in agreement with experimental NMR measurements. The confinement of the reacting molecules have also been taken into account by computations (Chakraborty and Chattaraj, 2018). However, these binding effects amounting to 1-2 kcal/mol, only partially explain the catalytic effects of 4-5 kcal/mol, which have so far been unresolved. We now investigate the reaction mechanism determining detailed underlying electronic structures to assess catalysis. Using continuum electrostatic model to account for the water environment, we obtained a good agreement between experimental and calculated reaction barriers without CB. However, in CB, the catalytic rate enhancement is not reproduced without taking into account explicit water molecules. Our QM/MM calculations reconcile these differences, and we are able to assess the catalytic effects upon binding in CB. Our results show that the catalytic effects are arising as a combination of confinement, electrostatic shielding and shifting the pK a of the ligand. In view of future design, we also demonstrate that DFT calculations can provide detailed electrostatic maps that elicit an optimal environment for efficient catalysis.

MATERIALS AND METHODS
Classical MD simulation were carried out using CHARMM36 force field (Best et al., 2012) as implemented in NAMD 2.12 (Phillips et al., 2005). A box of 3855 TIP3 waters solvating the complex was minimized in 20,000 steps, followed by 0.1 ns of NVT equilibration constraining the heavy atoms of 1a and CB [7]. Independent replicas of 1 ns unbiased MD were run with 2 fs timestep in NTP ensemble, coupled to a Langevin thermostat at a temperature of 303.15 K and a Langevin barostat at standard pressure. Non-bonded interactions were modeled with scaled 1-4 setting in NAMD and a cutoff of 12 Å with switching distance of 10 Å. The Particle Mesh Ewald method was used to evaluate long range electrostatic interactions. Subsequently, conformations for DFT calculation were chosen by conformation clustering done by a standalone utility shipped with Schrödinger Suite.
Intermediates and transition states have been optimized with B3LYP hybrid functional (Becke, 1993) and 6-31G * Pople basis set (Ditchfield et al., 1971;Hehre et al., 1972) using Grimme's D3 dispersion correction (Grimme, 2012) as implemented in Gaussian 09 Revision E (Frisch et al., 2009). Vibrational analysis and PCM solvent correction with Thrular's SMD parametrization was also carried out at the aforementioned level of theory (Marenich et al., 2009). For the sake of accuracy, final single point energies were calculated with 6-311++G(3df,3pd) basis set. Numerical integrals were computed on ultrafine grid. The nature of transition states was verified by Intrinsic Reaction Coordinate (IRC) calculations.
QM/MM calculations were done to study explicit solvation using TIP3P solvation model. The CB[7] and 1a was parametrized by CGenFF (Vanommeslaeghe and MacKerell, 2012), the guest was treated at B3LYP/6-31+G * level of theory. Instead of periodic boundary conditions, the system was trimmed to a sphere of 25 Å from the substrate, resulting in 2,907 or 2,866 explicit water molecules without or with CB[7], respectively. Atoms further than 20 Å away from the QM region were kept frozen throughout the simulations. Minimizations were run by CHARMM (Brooks et al., 2009) interfacing with Q-Chem 4.3 (Shao et al., 2015). The reaction path was sampled at discrete steps 0.2 Å apart along the reaction coordinate defined by the sum of the distances of the forming two sigma bonds. Quadratic constrains were defined with a force constant of 1,000 kcal mol −1 Å −2 , path scans were run until the convergence of the potential energy surface, then energies were recalculated using the 6-311++G(3df,3pd) basis set.
The cubic grid for the point charge analysis was made with a 0.5 Å spacing around the aligned reactant and transition state of the substrate 1a, to a maximum distance of 3 Å from the molecule. Points closer than 1.5 Å to any atom were removed. Barriers were recalculated with B3LYP/6-31+G * using Q-Chem 4.3.

RESULTS
We considered the cycloaddition reaction in four substrates (1a-d) obtained from substitution from N-allyl-2-furfurylamine molecule (Figure 1). To assess the catalytic effects, we compared reaction rates between water and in the presence of one equivalent of CB[7] (Figure 1B), which presents a 1:1 binary complex. To obtain the barriers, we identified reactant, product and transition states (RS, PS, and TS). We did not investigate the binding affinity of the substrates in CB[7], since the complexation is experimentally confirmed. As experimental reference results for the catalytic effects, we used barriers derived based on first order kinetics presented in reference (Palma et al., 2017). The introduction of CB[7] induces a ca. 5 kcal/mol drop in the activation free energy ( Table 1).

Protonation of the Substrates
First, we considered the effect of protonation on the amine moiety on the rate of cycloaddition. Since positively charged compounds are generally the most potent to bind into CB structures, due to the favorable interaction with the carbonyl  In experiments, pH was ∼7.4, while the pK a of the amine is estimated ∼10 based on the DFT free energies and empirical prediction (details in Supplementary Material). Henceforth, we consider a protonated substrate for all calculations. Nevertheless, we compared the neutral 1a compound in our preliminary calculations obtaining that it undergoes the cycloaddition with a barrier of 30.6 kcal/mol, while after protonation decreases this to 28.6 kcal/mol. This points to the beneficial effects due to the increased pK a of the amine in the CB environment, but does not fully account for the catalytic effects observed.

Reaction Without CB[7]
We calculated the reactant and transition states for all four differently substituted substrates to obtain the reaction barriers in solution. Geometry optimizations were carried out in the gas phase for the TS and for several reactant (Figure 2) and product conformers, followed by a continuum dielectric model to account for free energies in the water environment.
Considering the binding poses within the CB[7] cavity, the substrate needs to adopt a folded structure, which is readily aligned for cycloaddition (HIC, hair-pin induced conformation) FIGURE 2 | Conformers of substrate 1a without CB[7]: (A) most stable unfolded structure; (B) hair-pin induced (HIC) conformer directly preceding the Diels-Alder transition state (yielded by an IRC calculation). Relative Gibbs free energies are shown in kcal/mol. Non-polar hydrogens are hidden for clarity. (Palma et al., 2017). Therefore, we calculated the most stable open and HIC conformers of the reactant state for 1a that show a 1 kcal/mol free energy difference only. This demonstrates that additional effects beyond a HIC starting geometry also stabilize the transition state in the presence of CB.
Our calculated free energy barriers, summarized in Table 2, show a very good agreement with the experimental kinetic rates for all four compounds. The reaction barrier with the Me substituted 1b substrate is somewhat higher in the calculations as compared with experiment, whereas the other three ligands match the barriers within 0.2 kcal/mol. We also compared the free energy barriers using a range of functionals that all performed similarly. Small improvements can be observed using Minnesota functionals when we compare the barrier heights for the H and Me substituted ligands, 1a vs. 1b (Supplementary Table 2).

Reaction in the CB[7] Cavity
Initial MD simulations were used to obtain several conformers of host-guest complexes as starting points in DFT optimization for the stationary reactant, transition, and product states. Most We initially investigated the effects of a single water molecule for the catalytic barrier (Figures 3A,B). We placed the water molecule near the ammonium cation moiety given their strong interactions, and reoptimized the stationary structures. This resulted in small rearrangements in the binding mode of the substrate in the cavity, as the NH 2 does not directly interact with the CB[7]'s carbonyl oxygens, but through the explicit water instead (Figure 3C). The free energy barrier of the reaction in this model is 24.6 kcal/mol, which is in satisfactory agreement with the experimental results.
In comparison, analogous geometry optimizations for the non-catalytic system without CB[7], in the presence of a single water molecule does not change either the stationary structures or the barriers significantly (Supplementary Tables 1, 2). This suggests a cooperative influence of the solvent and the CB[7], which can only be incorporated with explicit solvation. Therefore, the stationary points along the reaction were calculated for all four substrates in the presence of one water molecule besides the macrocycle. The results are summarized in Table 2 showing a good agreement between this model and the  experiments. The level of theory for this series was benchmarked to cover other functionals and Pople basis sets, maintaining the same trends in all cases (Supplementary Tables 1, 2).

QM/MM Calculations in Explicit Water
To study the reaction featuring multiple water molecules, the most stable conformers of 1a were solvated in explicit water and subjected to reaction path scans using QM/MM calculations. The activation energy was 25.9 kcal/mol in solution, and 20.8 kcal/mol in complex with CB[7], minimized at the level B3LYP/6-31+G * and the energies recalculated with the basis set 6-311++G(3df,3pd).
The energy profiles depicted in Figure 4 exhibit a ca. 5 kcal/mol decrease in activation energy when the substrate is bound to CB[7], accounting for the overall catalytic effect fully. The only main difference in the reaction paths is observed for the reactant state regions. The reactant minima are located at reaction coordinate values 6.4 Å in water and 6.0 Å with CB [7]. Therefore, the RS conformations are more contracted inside the cavity before the cycloaddition, underlining the importance of pre-selecting conformers which easily undergo the cyclization. The overall barriers are somewhat lower than the ones presented in Table 2, which can be attributed to the different approach of solvent modeling.
Similarly, we also performed QM/MM calculations on the 1b substrate with a methyl substituent, as well as a nitro group (1e, Supplementary Figure 4) to analyze the trends in the barrier heights and the catalytic effects (Supplementary Figures 3, 5, respectively). We obtained a barrier height of 26.5 kcal/mol in water and 20.5 kcal/mol in CB[7] for 1b, showing a somewhat slower reaction in water, but faster in CB [7]. The QM/MM calculations thus improved the agreement with the experimentally observed catalytic effect for 1b, as also seen for compound 1a.
Overall, a clear trend is observed both in experiments and in calculations in water, without CB[7]. More electron withdrawing groups reduce the barrier more, therefore we predict the fastest reaction in water with the nitro group. Figure 5 interestingly, this trend is absent from the reactions with CB [7]. In this case, experimental barriers are comparable, no clear trend is observed.  Frontiers in Chemistry | www.frontiersin.org This is not observed in implicit solvent calculations with CB[7], and the prediction for the methyl substituent is affected the most. Therefore, the lack of explicit solvent in this case probably cannot accurately account for the ligand's geometry and bound position within the CB [7], which leads to the less accurate estimate of the catalytic effects. On the other hand, when using QM/MM, the presence of the water molecules allows us to obtain a significantly closer agreement with experiment for the catalytic effects. In this case, although the barriers both in water and in CB[7]-catalyzed environment are slightly underestimated, the overall catalytic effect is more accurately reproduced.

Point Charge Analysis
The well-defined Diels-Alder transition states enable us to study the shielding effect of CB [7] in an environment where molecular details are removed, and only electrostatic interactions are considered more generally. We placed probe point charges around the continuum solvent optimized reactant and transition states on a grid around the molecule and recalculated the energy barriers. We defined a cubic 3D grid up to 3.0 Å around the substrate with a 0.5 Å spacing to evaluate the electrostatic influence on the reaction barrier. Placing a −1 probe charge (given the cationic substrate) around the reacting molecule outlines the pattern on the catalytic effect depicted in Figure 6.
The surface was obtained by selecting the grid points that are located at within 2.0-2.3 Å distance from the ligand, where solvent molecules would surround the substrate. Yellow regions indicate an increased barrier height for the reaction when negatively charged groups, such as the oxygen of water molecules would be present, whereas green regions favor electrophilic groups at these locations for reducing the barrier heights. Our results using a +1 probe charge on the grid points results in analogous complementary effects (Supplementary Figure 6).

DISCUSSION
DFT calculations using continuum solvent models reproduce catalytic barriers to excellent accuracy as compared with experimental values. The results are consistent between tested functionals and basis sets within the acceptable error margin of DFT. However, the reaction barriers in binary complex with CB [7] are not able to capture the catalytic effects observed experimentally. Therefore, we used explicit water molecules as well as QM/MM calculations that successfully resolve the differences. Nevertheless, the absolute barriers are underestimated by the QM/MM minimizations, which can possibly be resolved by free energy calculation methods such as finite temperature string simulations. Our calculations suggest that the reaction path is slightly modified in CB [7] as there the compressed reactant state conformation is stabilized upon complexation. Explicit water molecules in this case enable us to obtain the correct binding orientation within CB [7], which lead to the decreased reaction barriers. Furthermore, electrostatic effects are also beneficial that shield the substrate from interactions with water.
Our results in conjunction with the experimental rates suggest that the four substituents have a similar catalytic effect (Figure 7). This is also seen from the similar geometric distances: C-C inner and C-C outer bonds observed in the stationary states determining the catalytic barrier heights. Therefore, these results can be generally applied to a range of molecules with analogous intramolecular Diels-Alder reaction mechanism that also contain an amine group, facilitating the binding and directing the orientation of the substrate. Extensions might include intermolecular reactions (Jon et al., 2001;Ren et al., 2018), where a larger CB[8] group might be necessary to accommodate both reacting partners.
The increasing electron withdrawing effects of the substrate substituents lead to faster reaction rates in water, as observed both in calculations and in experiment, including our predicted nitro substituent ( Figure 5). Interestingly, the catalytic effects also depend on geometric factors of these substituents, as the methyl and nitro groups deviate from the trend observed for the H, Br, and Cl series. The methyl group has a slightly enhanced catalytic effect, whereas the nitro group is relatively slower in CB [7].
To more generally assess electrostatic factors governing the reaction, we probed the catalytic effects of single point charges placed on a 3D grid around the substrate. We found that the regions, where the nearby negative charges unfavorably affect the reaction, largely overlap with the area that the CB[7] shields from solvent upon complexation. These results are consistent upon using positive probe charges that show complementary effects. We envision that these maps quantifying electrostatic effects for the reaction barrier will also enable future design for improving catalysis.

CONCLUSIONS
One of the first catalytic reactions where CBs were successfully used to enhance the rates of are the examples of Diels-Alder reactions. Our computational study elaborates the catalytic effect of cucurbit[7]uril on the cycloaddition reaction. The supramolecular host binds the substrate creating a favorably anisotropic environment preparing the cycloaddition both conformationally and electrostatically.
With DFT calculations and implicit solvation, we are able to reproduce the kinetics of the non-catalytic cycloaddition with great accuracy. However, the supramolecular catalytic effect is only interpretable by describing the environment more precisely, taking the solvent and their interaction with the substrate and the catalyst into account at the atomic level. Our results hence underline the limitations of implicit solvation and the cooperativity between the catalyst and the solvent.
In particular, a clear trend was observed in water for different substituents. The more electron withdrawing substituent used, the lower the barrier. This trend was not observed in the presence of CB[7] experimentally. In this case, the geometric effects due to interactions with solvent molecules likely play an additional important role in determining the reaction barriers, due to the ligand occupying the CB[7] cavity in slightly different positions. This is also supported by the fact that implicit solvent calculations continue to show a similar trend even in the presence of CB[7], whereas QM/MM calculations with explicit water, and computational models with added explicit water molecules more accurately reproduce the barrier heights for several substrates. Overall, both experimental and calculations suggest that the largest catalytic efficiency is achieved for the 1a and 1b substrates.
We recognize the natural direction of expanding the scope of this type of supramolecular catalysis based on the utilization the hydrophobic cavity and the H-bonding network formed by the solvent, polar groups and the portal of the CBs. Our general approach of mapping the influence of the electrostatics of the ligand's environment using simple 3D grid of point charges can aid this design of catalysts or additional substrates, where specific host-guest complex-based catalysis can expand to.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
OS and ER conceptualized the research, IS did the MD simulations, DB carried out the QM and QM/MM simulations and processed the results. All authors contributed to the article and approved the submitted version. FUNDING DB and ER acknowledge funding from ERC Starting Grant (Project 757850 BioNet). ER thanks funding from EPSRC (EP/R013012/1).

ACKNOWLEDGMENTS
IS and ER thank Ester di Tomasso and Benjamin Babiuck for their contribution. We are grateful to the UK Materials and Molecular Modeling Hub for computational resources, which is partially funded by EPSRC (EP/P020194/1). The authors acknowledge use of the research computing facility at King's College London, Rosalind (https://rosalind.kcl.ac.uk).