A DFT investigation of the host–guest interactions between boron-based aromatic systems and β-cyclodextrin

Density functional theory calculations including dispersion at BLYP-D3(BJ)/def2-SVP level of theory were performed for a series of systems based on cyclodextrin complexation with boron-based aromatic compounds. Elaborated investigations were carried out using different quantum chemical parameters such as computed complexation energies, theoretical association constants, and natural bond orbital (NBO) analysis. Several configurations and inclusion modes were considered in this work. The calculated complexation energies were consistent with the experimental classification of these systems on the basis of occurring interactions. Reduced density gradient (RDG) and independent gradient model (IGM) approaches determined the nature and strength of non-covalent interactions which played a central role in the formation of the complexes. Thus, phenylboronic acid (PBA) and benzoxaborole (Bxb) act mainly as hydrogen-bonded complexes with β-cyclodextrin, while mainly Van der Waals (vdW) interactions stabilize both catechol (PhBcat) and pinacol esters of phenylboronic acid (PhBpin) complexes. The ferroceneboronic acid (FcBA) exhibits a mixture of H-bonds and vdW interactions with β-cyclodextrin.


Introduction
In recent years, organoboron compounds have become one of the most versatile classes of heteroatom-substituted organic molecules [1][2][3]. Organoboron-based compounds and their derivatives are a class of organic molecules with several applications in various fields of chemistry such as analytical and supramolecular chemistry, catalysis, energy storage, and carbon capture (H 2 , CH 4 , CO 2 …) [4][5][6][7][8]. Indeed, boronic acids have been successfully used as catalysts, e.g., for the formation of amide bonds and enantioselective reactions [9][10][11][12]. Since 2005, boronic acids have been used by Omar Yaghi and his collaborators as secondary building blocks to synthesize nanoporous materials, in particular covalent organic framework like COF-1 [6].
The complexes which are based on the interaction of boric acid and its derivatives [13] with macromolecular systems such as carbohydrates have several important applications, notably in the area of drug delivery. For instance, it is reported that ferroceneboronic acid (FcBA) and its derivatives could act as electrochemical biosensors for detecting sugars [14].
Boron atom is more electropositive than carbon, and this fundamental property is exploited to the fullest in organic synthesis, becoming one of the most prominent areas of application of organoboron compounds. The mildly Lewis acidic character of the boron atom is suitable in the context of carbohydrate and fluoride sensing that is enabled by the formation of tetracoordinate borates with fluoride anions and polyols [3,[15][16][17][18][19].
In addition to the application areas mentioned before, boronic acids are used as therapeutic agents and biological probes. The first-in-class anticancer drug bortezomib has been joined by the recently approved antifungal tavaborole, and antieczema drug crisaborole, while many other boronic acid derivatives have shown promising activities in a large number of clinical and preclinical studies [20,21]. In this arena, boronic acids are complexed with supramolecular structures involving inclusion complexes with cyclodextrins. These studies have resulted in new commercially available drugs which have been of great benefit [22].
Cyclodextrins (CDs) have been used extensively due to their advantages of non-toxicity, facile modification, good water solubility, and high biological availability. CDs have been approved as pharmaceutical excipients for the manufacture of pharmaceutical preparations. CDs are cyclic oligosaccharides formed by 6, 7, or 8 glucose units, which are called α-, β-, and γ-cyclodextrins, respectively (see Fig. 1) [23][24][25].
The 3D structure of CDs shows hollow truncated coneshaped macrocycles in space; their inside cavity is hydrophobic while the outside is hydrophilic. CDs form inclusion complexes with a wide variety of boronic acids [26][27][28][29]. In this context, Kasprzak et al. studied the interactions of some boric acid derivatives with β-cyclodextrin (β-CD) using spectroscopic methods: proton nuclear magnetic resonance ( 1 H NMR), 1H diffusion-ordered spectroscopy nuclear magnetic resonance ( 1 H DOSY NMR), and Fourier-transform infrared spectroscopy (FTIR) [30]. They found that phenylboronic acid (PBA) and benzoxaborole (Bxb) form non-covalent hydrogen bondingbased systems with β-cyclodextrin, whereas catechol (PhBcat) and pinacol esters of phenylboronic acid (PhBpin) as well as ferroceneboronic acid (FcBA) form host-guest inclusion complexes. It has been found that the guest molecules Bxb and PBA which contain free hydroxyl groups are likely to form noncovalent systems based on the hydrogen bonding with β-CD, due to the cross-correlations observed in the 1 H-1 H ROESY NMR spectrum between O-H groups of Bxb and PBA and the O-H groups of β-CD; however, the guest molecules PhBcat and PhBpin are moving inside the β-CD inner cavity through their aromatic rings, and cross-correlations between H-3 and H-5 of β-CD and aromatic protons of PhBcat and PhBpin were observed in the 1 H-1 H ROESY NMR spectrum, thus proving their complexation inside β-CD cavity and leading, therefore, to the formation of inclusion complexes. The proposed structural models for the studied complexes are represented in Fig. 2.
To quantify the interaction between the studied boronbased aromatic systems and β-cyclodextrin, Kasprzak  Furthermore, despite the diversity of characterization techniques used, the authors did not explain the mechanism of inclusion. That is why the current study is aimed at a more profound understanding of the inclusion complex formation between select aromatic boron-based compounds (PBA, PhBcat, PhBpin, Bxb, and FcBA) and β-CD using density functional theory. To the best of our knowledge, this is the first theoretical study aiming to quantify the interactions between a class of aromatic boron compounds and β-cyclodextrin using exclusively DFT-based calculations.

Computational details
DFT-based computations were carried out with ORCA program (version 4.2.0) [31] and Gaussian 09 quantum computational package [32]. The geometry of each resulting configuration was fully optimized in both gas phase and aqueous solution (CPCM solvent model) at BLYP-D3/def2-SVP level of theory [33] in conjunction with Becke-Johnson (BJ) damping function dispersion correction D3(BJ) [34][35][36][37][38] and def2-SVP basis set [39][40][41]. For improving the computational efficiency of our large-scale calculations, we employed the resolution of the identity (RI) approximation [42,43]. The geometries of β-CD, guest molecules, and the most stable formed complexes were optimized and verified to be stationary points by frequency analysis. A geometrical counterpoise correction scheme (gCP) to the def2-SVP for the intra-and inter-molecular basis set superposition error was applied [44]. Default values convergence criteria were used for all optimizations. Several studies have shown the reliability of BLYP-D3(BJ) in describing satisfactory non-covalent interactions [45][46][47][48].
The natural bond orbital (NBO) calculations were carried out with the Gaussian 09 code in the same manner as the precedent protocol (BLYP-D3(BJ)) but with the larger basis set def2-TZVPP to elucidate the nature of hydrogen bonding in the studied complexes.
The initial structure of β-CD was taken from Chem-Office 3D ultra (Version 10, Cambridge Software) [49]. The starting geometries of studied boron compounds were constructed using Hyperchem 7.5 molecular modeling package [50].
The coordinate system defining the process of complexation is represented in Fig. 3.
Following the method proposed by Liu and Guo [51], the glycosidic oxygen atoms of β-CD were positioned onto the XY plane; their center was set as the origin of the coordinate system. The secondary hydroxyl groups of the β-CD were placed pointing toward the positive Z-axis. The guest molecules (aromatic boron compounds) were initially placed at the center of the coordinate system (0 Å) and then translated from −5 to +5 Å along the Z-axis in both directions with steps of 1 Å, thus resulting in two possible modes of inclusion A and B (Fig. 3). The structures of guest molecules were initially pre-optimized. All the generated configurations were fully optimized in both gas phase and aqueous solution (CPCM solvent model) at BLYP-D3/def2-SVP level of theory. The complexation energies for each structure were determined using the following equation: where ΔE complexation represent the energy gain due to complexation and E guest@−CD ,E guest@ −CD and E −CD are the energies of optimized geometries of the complex, the free guest, and the free β-CD, respectively.
For an exploration of the nature and strength of the interactions [52][53][54][55][56] existing between the aromatic boron compounds and β-CD, we performed a reduced density gradient (RDG) [57] and independent gradient model (IGM) analysis [58]. The non-covalent interaction (NCI) [59] maps were characterized with Multiwfn software [60] and visualized with VMD program [61]; they are represented by generating colored graphs of RDG isosurfaces, where the blue, green, and red regions are associated, respectively, to H-bonds, Van der Waals interactions, and steric effect.

Inclusion complexation energy evaluation
The values of the computed complexation energy in gas and aqueous phases as a function of the Z coordinate during the inclusion process for A and B models are reported in Tables 1 Fig. 2 Interactions of boronbased aromatic systems with β-cyclodextrin [30] and 2, where both models exhibit negative energies indicating that the occurring process is thermodynamically favored.
When the guest molecules approach β-CD in gas-phase calculations, several local minima were observed, but the most sta-  (Table 1).
The complexation energies are more negative in vacuum than in aqueous solution; they range between -147.4 and -172.2 kJ/mol in gas phase and between -83.8 and -105.2 kJ/mol in water solution. Due to the water solvent effects, the interaction between aromatic boron compounds and β-CD is weakened; showing, therefore, that inclusion process is more exothermic in gas phase than in water.
The calculations carried out in aqueous phase show the same trend toward the complexation energy except that FcBA@β-CD is more stable than PhBpin@β-CD.
The dispersion interactions can play a key role in stabilizing non-covalent complexes. In order to analyze their contributions to the total interaction energies for the studied complexes, single-point calculations at M06-2X-D3(0)/TZVP level [62] were performed on geometries optimized at BLYP-D3BJ/def2-SVP. It is well established that M06-2X is one of the best functionals for the study of non-covalent interactions [45,46,63,64].
The Overall, the dispersion energies of the studied complexes are correlated with the experimental nature of the assembly between the boron-based aromatic systems and β-CD. Indeed, the PhBcat@β-CD, PhBpin@β-CD, and FcBA@β-CD can form inclusion complexes, whereas PBA@β-CD and Bxb@β-CD favor the formation of hydrogen-bonded systems.
The resulting gas-phase optimized geometries at BLYP-D3(BJ)/def2-SVP level will be used in the subsequent calculations of NBO and RDG function.

Theoretical determination of association constant
From a computational viewpoint, the association constant K a could be calculated according to the following reaction: Fig. 3 Coordinate system used to define the inclusion process between studied boron compounds and β-CD. Atom's color code: pink for B, purple for Fe, grey for C, red for O, and white for H However, the association constants are determined experimentally in the aqueous phase. To be able to compare the experiment K a with the results of quantum chemistry calculations, it is, therefore, necessary to take into account the effects of solvation on the calculation of the free reaction enthalpies using CPCM solvent model.  It is known that in these types of systems, a continuous model for the solvent can strongly fail since it does not allow evaluating the role of specific interactions between the water molecules and the guest/host molecules. Nevertheless, Champion et al. have proposed a strategy and applied it successfully in predicting reaction equilibrium constants, for instance, for reactions similar to the ones studied in our work [65,66]. They have shown that more accurate results can be obtained when the thermodynamic cycle involves not direct complexation reactions, but rather ligand-exchange reactions by determining the exchange constant Kexc. The basic reasoning justifying the use of this strategy is based on error cancelation mechanisms between the species on the left-and the right-hand side. We hence proceeded accordingly, and, in our calculations, we employed boron compounds exchange reactions involving the guest@βCD species as shown in Fig. 4. We chose the phenylboronic acid catechol (PhBCat) as a reference.
The boron compounds exchange equilibrium constant, K exc , is then calculated from the free energy change in aqueous solution, ΔG * Sol , which can be expressed as: where G 0 g is the free energy change in the gas phase and ΔG * sol (guest@βCD) , ΔG * sol (PhBcat@βCD) , ΔG * sol (guest) , and ΔG * sol (PhBcat) are the free energies of solvation of the respective species in water. For each species, its Gibbs free energy at 25 °C is obtained using the computed energy with the BLYP-D3(BJ)/def2-SVP level, and the thermodynamic corrections from the frequency calculation were performed with the same functional and basis set in gas phase. The computed Gibbs free energies at 25 °C are regrouped in Table S1. Note that some guests and guest@β-CD complexes exhibit several competitive configurations (see, for instance, the case of Bxb@β-CD; Fig. S1); their Gibbs free energies have then been evaluated using a Boltzmann distribution: where the summation runs over all the most stable configurations of the A species. The cartesian coordinates of the different configurations used in the determination of the exchange constant are presented in Table S2. The exchange equilibrium constants (Log K exc ) computed at the BLYP-D3/def2-SVP level of theory, and the experimental values, are presented in Fig. 5.
From the computed Log K exc (the dots in green color; Fig. 5), it is clear that the complexes PhBcat@β-CD, PhBpin@β-CD, and FcBA@β-CD (inclusion complexes) have the three strongest association constants, while the complexes forming hydrogenbonded systems have the lowest ones. This finding confirms the previous results concerning the complexation energies.
On the other hand, we notice that for the FcBA, PhBcat, PhBpin, PBA, and Bxb, the computational values agree well with the experimental results. The mean absolute deviation (MAD) between the experimental and computed Log K exc values is rather small (0.37).

Intermolecular hydrogen-bonding effects
This section is intended to shed light on the ability of the aromatic boron compounds to interact through hydrogen bonds (HBs) with β-CD.
NBO method consists of interpreting the electronic wave function in terms of Lewis structures by considering all possible interactions between filled donor (i) and empty acceptor (j) NBOs and evaluating their stabilizing energy through the second-order perturbation theory. The results of NBO analysis associated with the bond length and the stabilization energy of intermolecular hydrogen bonding in the studied complexes are reported in Table 3. Fig. 4 Thermodynamic cycle used to compute free energy changes in aqueous solution  A comparison of significant hydrogen bond lengths shows that the values fluctuate between 1.85 and 2.76 Å (Fig. 6), and the shortest lengths (1.85-2.10 Å) corresponding to strong donor-acceptor interactions [67,68] with higher stabilization energies are associated with the complexes PBA@β-CD, Bxb@β-CD, and FcBA@β-CD.
The strength of intermolecular hydrogen bonds varies from weak for PhBcat@β-CD and PhBpin@β-CD to moderate for the complexes PBA@β-CD, Bxb@β-CD, and FcBA@β-CD. Indeed, as reported in Table 3, the intermolecular hydrogen bond lengths of PhBcat@β-CD and PhBpin@β-CD complexes are averagely longer than those of PBA@β-CD, Bxb@β-CD and FcBA@β-CD and their corresponding energies are consequently lower.

Nature of non-covalent interactions
RDG and IGM plots were carried out by an isosurface value of 0.50 and illustrated in Fig. S2 (supplementary data). Red, green, and blue colors represent respectively steric repulsion, weak Van der Waals interactions, and strong attractive interactions like H-bonds.
The RDG results illustrated in Fig. S2 show clearly that the inclusion of PBA, PhBcat, PhBpin, Bxb, and FcBA molecules in β-cyclodextrin exhibited different behavior. The two guest molecules PhBcat and PhBpin are wrapped by green isosurfaces, indicating weak Van der Waals interactions occurring between the PhBcat Visual inspection of the RDG scatter plots shows the new spikes compared with isolated β-CD (Fig. 7), in the ranges [0.000 to -0.015] and [-0.020 to -0.025] a.u. which correspond, respectively, to the vdW and hydrogen bond intermolecular interactions. Additionally, the intermolecular interactions between β-CD and guest molecules are quantified through the independent gradient model (IGM) to identify the nature of the intermolecular interactions, which can be provided by the δginter and δgintra. Figure 8 represents the IGM isosurfaces of δginter for the five complexes, where the green-colored isosurfaces denote weak van der Waals interactions, whereas blue regions indicate stronger attractive interactions associated with hydrogen bonds.
The scatter plots between δginter and δgintra versus sign(λ2)ρ were shown in lower figures of Fig. 8, in which the red and black points correspond respectively to δg inter and δg intra .
The plotted IGM isosurfaces confirmed the nature of occurring interactions. Indeed, the complexes PBA@β-CD and Bxb@β-CD can be considered as hydrogen-bonded complexes, while mainly vdW interactions stabilize PhBcat@β-CD and PhBpin@β-CD. The formation of the complex FcBA@β-CD is a combination of H-bonds and vdW interactions.
Overall, van der Waals (vdW) dispersion interactions can strongly affect the stability of the formed complexes.

Conclusions
In conclusion, we proposed the first DFT study of five systems based on cyclodextrin complexation with a series of aromatic boron compounds. The DFT-D3 investigation results confirmed their experimental ranking based on the formation of inclusion complexes or hydrogen-bonded systems and showed that dispersion interactions are the major contributors to the formation of complexes.
The inclusion process performed in vacuum and water showed that PhBpin@β-CD and FcBA@β-CD have respectively the highest complexation energies.
The natural bond orbital (NBO), reduced density gradient (RDG), and the independent gradient model (IGM) calculations in gas phase showed that PBA@β-CD and Bxb@β-CD form hydrogen-bonded systems, while mainly vdW interactions stabilize PhBcat@β-CD and PhBpin@β-CD. FcBA@β-CD complex is stabilized by a combination of H-bonds and vdW interactions.
One significant result of this study is the higher stability of ferroceneboronic acid upon its inclusion in the β-cyclodextrin cavity; since ferroceneboronic acid and its derivatives are used as biosensors, this study could serve as a starting point in experiments for overcoming problems of solubility and bioavailability by considering their inclusion in β-cyclodextrin.