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., 2017, 2019; Olesiń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., 2017, 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, 2014; Lee 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.

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⁻¹ Å⁻², 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.

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

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.

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.

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

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.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.