Unveiling the gemcitabine drug complexation with cucurbit[n]urils (n = 6–8): a computational analysis

In this work, we have investigated the complex formation capability of gemcitabine drug with host cucurbit[n]urils, Q[n] (n = 6, 7, and 8) using density functional theory (DFT). The DFT studies demonstrate that the most stable configuration is a fully encapsulated complex. In the encapsulated gemcitabine@Q[6] and gemcitabine@Q[7] complexes, the gemcitabine amino -NH2 and the alcoholic groups bond with the carbonyl units of Q[n]. The addition of sodium ions leads to the partial exclusion of the gemcitabine molecule and the sodium atoms lie close to the carbonyl portal of Q[7]. Thermodynamic parameters computed for the complexation process exhibit high negative Gibbs free energy change, implying that the encapsulation process is spontaneous and is an enthalpy-driven process. Frontier molecular orbitals are located mainly on the gemcitabine uracil ring, before and after encapsulation, indicating that the encapsulation occurs by pure physical adsorption. Quantitative molecular electrostatic potentials demonstrate a shift in electron density during the complex formation and are more pronounced in gemcitabine@Q[7]. Atoms-in-molecule (AIM) topological analysis illustrate that these complexes are stabilized by various non-covalent interactions including hydrogen bonding (HBs) and C···F interactions. The reduced density gradient (RDG) graph exhibits the presence of strong HBs and weak van der Waals interactions and the presence of steric repulsion. The isosurface non-covalent interactions (NCI) diagram shows predominant steric interaction in the gemcitabine@Q[6] complex. The NCI isosurface for gemcitabine encapsulated complexes with Q[7] and Q[8] host displays that the green patches are uniformly distributed in all directions. Finally, EDA results demonstrate Pauling’s repulsive energy is predominant in gemcitabine@Q[6] complex, while the orbital and dispersion energies stabilize the gemcitabine@Q[7] complex.


Introduction
Macromolecules are repeatedly appealing and fascinating due to their interesting topology structures having repeated subunits with excellent host-guest properties [1,2]. Among the known macrocyclic molecules, cucurbiturils find a prominent position as their synthetic techniques lead to molecules that are soluble in aqueous as well as in organic solvents [3]. Furthermore, it was found that Q[n] (n = 5-9, 11, and 14) binding with both inorganic and organic molecules. The binding constant of Q[n] with ferrocene was up to 10 15 M −1 , nearly equivalent to that biotin-avidin pair which is considered to be the system with the highest binding biomolecular system identified by humans [4]. These synthetic receptors are found to have a simple structure with a pumpkin-like shape. The repeating units found in these molecules are the cyclic methylene groups bridged with the glycoluril oligomers with two portals on the top and bottom by the ureido carbonyl groups [5]. The formation of the cyclic portal with a carbonyl group at the top and bottom leads to the formation of a hydrophobic cavity in the middle of the molecule, offering a unique property for the attraction of guest molecules with high binding efficiency.
Gemcitabine (2′,2′-difluorodeoxycytidine) is one of the most widely used in cancer treatment and third in clinical oncology [6,7]. In specific, it is mainly used in pancreatic adenocarcinoma treatment and is also prescribed to treat a variety of solid tumors including breast, ovarian, and nonsmall cell lung cancers [8]. Gemcitabine has a structural similarity with deoxycytidine, with the former having two fluorine atoms replacing the hydrogen atom at the C2′ carbon. Gemcitabine when administered in the human body is either deaminated by deoxycytidine deaminase in the liver to inactive 2′,2′-diflurodeoxyuridine or activated by deoxycytidine kinase to 2′,2′-diflurodeoxyuridine-5′-monophosphate [9,10]. Furthermore, 10% of the gemcitabine undergoes renal filtration, and within a week, more than 90% of the dose gets removed as 2′,2′-diflurodeoxyuridine in the urine [11]. Due to the hydrophilic nature of the glucose ring, gemcitabine cannot infiltrate into cells by passive diffusion [12]. Previous clinical studies indicate that lower doses can be effective and it could be advantageous to deliver gemcitabine [13]. Consequently, the delivery of gemcitabine in a lower dose, with prolonged systematic exposure with flexibility in administration using an oral formulation, is a desire for its wide use.
Cucurbiturils show comparatively low toxicity and make guest molecules crossing through the cellular membrane easier [14,15], hence are widely used as a drug delivery vehicle for both metallo and organic pharmaceuticals. In the past, cucurbiturils have been used as drug delivery vehicles for cisplatin, proflavine, and doxorubicin [16][17][18][19]. Buczkowski et al. have carried out an in vitro cell viability study on tumor cell lines MOLT4, THP-1, and U937 and found that free gemcitabine and gemcitabine complexed by Q [7] show equal cytotoxicity when exposed for 24 h and 72 h [20]. A 1:1 stoichiometry was observed at a low concentration of Q [7], while a molar excess of cucurbituril results in a 1:2 complex formation. Thermodynamic studies show the reaction was exothermic and is accompanied by a change in entropy [21]. Very recently Adam Buczkowski studied the binding of gemcitabine with Q [7] in the pH range of 5 to 2. The addition of NaCl decreases the binding rapidly, which was linked with the binding of sodium ions to the Q [7] carbonyl portals.
The present study intends to divulge the nature of interactions and thermodynamics of the complexation of three cucurbituril homologs (Q [6], Q [7], and Q [8]) with the gemcitabine drug in the gas and solution phase. A series of DFT calculations were employed to identify the stable conformers of the cucurbituril complexes with gemcitabine molecule. The role of Na ions in the complexation was found by introducing one and two Na atoms to the gemcitabine molecule. To address the queries concerning the factors controlling the complexation, we computed the thermochemical parameters in aqueous media. To understand the nature of interactions responsible for the stability of the complexes, quantitative molecular electrostatic potential (MESP), atoms-in-molecules (AIM), noncovalent interaction (NCI) analysis, and energy decomposition analysis (EDA) are carried out. The outcome of the current study would assist in the encapsulation of new drug molecules and their use as drug delivery vehicles.

Computational details
The initial geometries of gemcitabine and cucurbit[n]urils Q[n] (n = 6-8) are adopted from previously published results [23,24]. The structures were optimized without any symmetrical constraints using M06-2X functional including Grimme's empirical dispersion D3 correction [25]. Previous benchmark studies have shown that M06-2X-D3 functional delivers the best agreement for their geometry and bond length with crystal structures [26]. Furthermore, Assaf et al. suggest that the use of DFT-D on the host-guest complexes and macromolecular complexes improves the binding energy of the complex [27]. We used the 6-311G(d,p) basis set in all the calculations. The use of 6-311G(d,p) basis set is vindicated as a compromise between acceptable results and reasonable computational cost [28]. The solvent effects were applied through the conductor-like polarizable continuum method (CPCM), with sphere radii optimized for COSMO-RS proposed by Klamt [29]. Vibrational frequencies were computed for all stated configurations to confirm them to be real minima on the potential energy surface and to compute the thermodynamic parameters. All the DFT calculations were carried out using the Gaussian G16 suite of program packages [30]. The binding of the gemcitabine molecule is calculated using the supramolecular method [31].
The wave function analysis-surface analysis suite (WFA-SAS) program was employed to compute the quantitative electrostatic potential at 0.001 electron/Bohr 2 isodensity and to picture the 3D surface [32]. The M062X-D3/6-311G(d,p)//M062X-D3/3-21G method was employed for the generation of wavefunction to compute the electrostatic potential surfaces. The EDA is carried out at the M062X-D3/6-311G(d,p)//B3LYP-D3/TZ2P level using the Amsterdam density functional (ADF) program package [33]. In EDA, the total binding energy (ΔE bind ) can be disintegrated into four terms, viz., Pauli repulsion (ΔE Pauli ), electrostatic (ΔV elstat ), orbital (ΔE orbit ), and dispersion (ΔE disp ) terms [34]. For performing EDA analysis, the gemcitabine is considered as one fragment and the Q[n] as the other fragment. The repulsive term ΔE Pauli accounts for the repulsive Pauli interaction between the occupied orbitals of two interacting fragments. The attractive terms ΔV elstat and ΔE orb represent the classical electrostatic interactions between the fragments and the interactions between the occupied molecular orbitals of one fragment with the unoccupied molecular orbitals. The term ΔE disp represents the dispersion interaction among two interacting fragments. AIM analyses were carried out using the AIMALL program package on the wavefunction obtained on M062X-D3/6-311G(d,p) method [35]. NCI analyses were accomplished using MultiWFN software and visualized using the CHEMCRAFT program [36,37].

Structure of host, guest, and complexes
The chemical structures of cucurbit[n]urils, Q[n] (n = 6, 7, and 8) and gemcitabine are shown in Fig. 1a, b respectively. The drug gemcitabine is a pyrimidine nucleoside prodrug in which the hydrogen atoms on the 2′ carbon of deoxycytidine are replaced with fluorine atoms. The computed bond parameters of Q[n] were found to agree well with the previous reports [38]. The most important parameters used to judge the ability of a host to transport a well-defined drug molecule through the boding are its binding energy along with its ability to encapsulate the drug. To identify the best conformer and to know the possible encapsulation of the drug inside the Q[n] (n = 6, 7, and 8), the drug is allowed to interact with the Q[n] in all possible modes. Geometry relaxation was carried out both in the gas and solution phases. In the solution phase calculations, we employed the conductor-like polarizable continuum model (CPCM) using the COSMO-RS radii. We noticed several adsorption configurations including the full, partial, and surface adsorption sites for the drug in all the Q[n] molecules. The partial encapsulation includes either encapsulation of the pyrimidine base or the carbohydrate part. The most stable geometry for the guest complexation is shown in Fig. 1c-h. The other low-lying gemcitabine complexes with Q [6], Q [7], and Q [8] are shown in supporting information Figs. S1-S3.
It is interesting to observe that the most stable configuration is encapsulated complexes both in the gas and solution phases. Furthermore, gemcitabine is fully encapsulated in Q [6] and Q [7], and the amino -NH 2 and the alcoholic group in the carbohydrate of gemcitabine bonds with the carbonyl units of Q [6] and Q [7]. However, the guest orientation differs by the size of Q[n]. In the gemcitabine@Q [6], the gemcitabine is perpendicular to the Q [6]. In the gemcitabine@Q [7] and gemcitabine@Q [8], we observe the guest orientation to be partially and fully parallel to the carbonyl portal of the hosts. Our calculations show that gemcitabine shows a packing coefficient (PC) of 82, 66, and 72% with Q [6], Q [7], and Q [8]. Among the hosts, Q [7] has the lowest PC value of 66% and the highest for Q [6]. This clearly shows that in Q [6] the gemcitabine undergoes a tight fit into a host, while Q [7] has a perfect packing coefficient needed for the stabilization of an inclusion complex as suggested by Rebeck Jr. [39]. However, recent studies by Grimm et al. have identified the hostrelated solvation effects to be the major factor controlling the affinity of guest molecules with Q[n] complexes [40].
During complexation, structural changes in the host and guest were observed and are found to be the least for the most stable complexes. The strain energy is computed from the difference in energy values of the host and guest in their ground state geometry and the energy obtained on the geometry after complex formation. The computed strain energy for the encapsulated complexes in the solution phase is shown in Table 1. The computed strain energy shows large variation in the guest as well as the host molecule in the encapsulated drug-host complexes. The strain energy on the host and guest molecules was the minimum for the gemcitabine@Q [7] and highest for the gemcitabine@Q [6]. Recently, Athare and Gejii suggested that for the most stable complexes, the deformation energy will be minimum on the host and guest molecules [41]. In the present study, we observe the structural changes were minimum in the host and drug molecule for the gemcitabine@Q [7] encapsulated complex, which supports the previous findings [42].
A close observation shows that the hydroxyl groups of the deoxycytidine and the amino groups are found to exist in HBs with the carbonyl unit of the Q[n] molecule. Interestingly, in the encapsulated complexes, the number of HBs varies with the complex. In the gemcitabine@Q [6] and gemcitabine@Q [7] complexes, two hydroxyl groups and one amino hydrogen of gemcitabine were bonding with the cucurbiturils carbonyl group. In the case of gemcitabine@Q [8], one hydroxyl group and two amino hydrogens were bonding with the cucurbiturils carbonyl unit. In the gemcitabine@Q [6] stable complex, the two hydroxyl groups are at a distance of 2.097 and 2.178 Å from the carbonyl group of Q [6]. In the most stable gemcitabine@Q [7] complex, the above distances are 2.163 and 1.942 Å, while in the gemcitabine@Q [8] the distance was 1.847 Å. Likewise, the amino groups were found to have different bonding distances in the complexes. Thus, in general, the least deformation of the host and guest molecules exemplifies the host as an optimal fit for the guest and the presence of HBs  [6], d top perspective view of gemcitabine@Q [7], e top perspective view of gemcitabine@Q [8], f lateral view of gemcitabine@Q [6], g lateral view of gemcitabine@Q [7], and h lateral view of gemcitabine@Q [8] obtained from M062-X/6-311G(d,p) method in water between gemcitabine and carbonyl portal of the Q [7] host makes gemcitabine@Q [7] a stable complex.
To understand the effect of pH on the gemcitabine drug complexation with Q [7], recently Adam Buczkowski studied their interaction in both acidic and aqueous solutions [21,22]. He observes that a decrease in pH improves the binding of the drug to Q [7] and an increase in sodium chloride concentration decreases the complex stability, presumably due to the hindering effect of the binding of sodium cations with Q [7] portals. To know the real trend, we introduced one and two sodium atoms by replacing the hydrogen atoms of the hydroxyl group as they are more prone to basic conditions. The fully optimized geometries of the inclusion complex are shown in supporting information Fig. S4. After the introduction of sodium, we observe that in an encapsulated complex, the amino group of gemcitabine is moved out of the Q [7] cavity. With the introduction of the second sodium, we observe that gemcitabine moves out of the cavity more. Furthermore, the formed complex loses the HBs between the gemcitabine and Q [7] and the sodium atoms lie close to the carbonyl portal of Q [7].

Energetic analysis and stability of complexes
To know the stability of the encapsulated gemcitabine with Q[n] (n = 6, 7, and 8), we calculated their binding energy. The binding energy is computed as the energy difference between the total energy of the gemcitabine encapsulated complex and the total energy of various cucurbituril and gemcitabine molecules at their most stable geometry [42]. The binding energy of gemcitabine with various cucurbiturils in gas and solution phases is provided in Table 1. The BE values are found to vary with the Q[n]. Furthermore, gas-phase BE values are found to be very high compared to the solution phase. The higher values in the gas phase calculation ascend as the water molecules in the cavity of Q[n] were not considered in the gas phase simulations (which are regarded as high-energy water). While in the solution phase, those water molecules make hydrogen bonding with guest and host molecules and stabilize the complex [43][44][45]. For the determination of accurate binding energies, highenergy water molecules must also be considered as they contribute directly to the cavitation energy [46,47]. Very recently, Grishaeva et al. have revealed that high-energy water molecules play a more decisive role in the anchoring of guest molecules inside the cavity of Q [8] than in Q [6] and Q [7] host [48]. Among the inclusion complexes in the solution phase, gemcitabine@Q [6] has the least BE value of −28.57 kcal mol −1 , while gemcitabine@Q [8] has the −39.87 kcal mol −1 . The largest BE value was observed for the gemcitabine@Q [7] with −42.14 kcal mol −1 . Thus, among the hosts studied, gemcitabine has the highest binding energy with Q [7] and the least with Q [6] molecule. The computed binding energies of gemcitabine with Q [7] in the presence of one and two sodium ions were found to be −41.95 and -43.15 kcal mol −1 , which are close to the gemcitabine@Q [7] complex. It is worth pointing out that in these calculations the implicit solvation alone is considered, while the presence of explicit solvation is essential to stabilize the cations [47]. If the explicit solvation is taken into account, then computed binding energies would be much lower than the reported values. In the ITC experiment, the equilibrium constant value was found to decrease with the increase in the sodium ion concentration which corroborates our findings [22].
To deduce the complexation capability and dependency of gemcitabine with hosts Q[n] (n = 6, 7, and 8), we calculated the thermodynamic parameters such as standard free energy change (ΔG°), enthalpy change (ΔH°), and entropic component (TΔS°) for these encapsulated complexes in solution phase. Table 1 summarizes the thermodynamic parameters for the encapsulated complexes calculated at 1 atm and 298.15 K using the partition functions obtained in a vibrational frequency calculation. Adam Buczkowski reported the ΔH°, ΔG°, and TΔS° for the complexation of gemcitabine with Q [7] from an isothermal titration calorimetry (ITC) experiment at pH = 4.5 as −41.8 ± 1.6, −20.6 ± 0.2, and −21.2 ± 1.7 kJ mol −1 respectively. Our reported values are an order less than the experimental values reported by Adam Buczkowski [21]. The enthalpy changes were negative for all the studied complexes, implying the encapsulation process is exothermic. Besides, the Gibbs free energy change was also negative, demonstrating the encapsulation process is spontaneous and thermodynamically favorable. The high negative ΔG° and ΔH° values suggest that the encapsulation process is spontaneous and is an enthalpy-driven process [49,50].
The computed entropies were also negative indicating that the system transfers to a more ordered state during encapsulation compared to the free fluctuating guest/host molecules. It is interesting to find that the ΔG° and ΔH° are higher for the Q [7] encapsulated complex than Q [6] and Q [8] complexes. On the contrary, the entropy contribution for all the complexes was more negative than other thermodynamic parameters and are nearly similar, implying that for the complexation process enthalpy plays a major role than entropy. To understand the effect of temperature during the complex formation process, we computed the ΔG° at temperatures 298, 273, and 216 K which are listed in Table 1. With the decrease in the temperature, we observe that the complex formation increases. Thus, from the thermodynamic point of view, the formation of gemcitabine@Q [7] encapsulation is more facile, stable, and spontaneous compared to the other studied complexes.
Previous theoretical studies correlate the HOMO-LUMO gap as an index of kinetic stability for the molecules [51]. To know the kinetic stability of these encapsulated complexes, we calculated the energy separation of the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). The HOMO-LUMOs gaps (E gap ) computed for the encapsulated complexes are listed in Table 2. The Q[n] molecules have an E gap of 9.75-9.75 eV, while the gemcitabine has E gap of 7.93 eV. The inclusion complex has E gap of 7.70-7.91 eV which is lower than gemcitabine and the host Q[n] molecules. This shows that the gemcitabine molecule is chemically stable when encapsulated and the encapsulation process is reversible and kinetically controllable.
The frontier molecular orbital (FMO) of the encapsulated gemcitabine with cucurbiturils is provided in Fig. 2, while the hosts and guest are provided in the supporting information in Fig. S5. From Figs. 2 and S5, it is evident that the HOMO and LUMO orbitals are located mainly on the gemcitabine uracil ring, before and after encapsulation formation, indicating that the encapsulation does not alter the electronic nature the of drug and the process occurs by pure physical adsorption [52].
From HOMO and LUMO eigenvalues, we computed global reactivity descriptors (GRD) such as electronic chemical potential (μ), chemical hardness (η), and global electrophilicity index (ω) which are connected with the maximum hardness and minimum electrophilicity principle (MHP). The above quantities are used to arrive at electrophilicity-based charge transfer (ECT). The details to attain GRD parameters are provided in our previous publications [24,53] and the values obtained are listed in Table 2. The chemical potentials of all the complexes were negative, which suggests that the encapsulated complexes are stable. Also, the chemical potentials of the gemcitabine molecule are greater than the free host molecules and the encapsulated complexes have chemical potential higher than the free guest and lower than the gemcitabine. The chemical potential indicates the measure of electron transfer in a system during the complexation process. The chemical hardness of the gemcitabine molecule is larger than the cucurbit[n]uril molecules. Among the cucurbit[n]uril molecules, no appreciable change in hardness values was observed. Also, the chemical hardness values of the encapsulated complexes are closer to the free guest molecule. The electrophilicity index ω is a positive number and larger values are characteristics of most electrophilic systems [54]. From the table, we can see the guest molecule has higher ω values than the host molecules. Furthermore, among the encapsulated complexes, gemcitabine@Q [7] has a lower electrophilicity index ω, suggesting that Q [7] is a better host for gemcitabine. We have calculated the electrophilicity-based charge transfer (ECT) using Eq. (1) where ΔN max = −μ/η. The gained results for the encapsulated complexes are listed in Table 2. The ECT values were positive for all the complexes, demonstrating the charge transfer occurs from the guest to the host molecule and that the charge transfer from the guest is nearly identical for all the host molecules. Thus, in addition to the charge transfer, other intermolecular interactions also show a vital role in stabilizing the gemcitabine inside the host molecules.

Nature of intermolecular interaction
To realize the charge transfer from the host to the guest and to identify the other intermolecular interactions that are accountable for the stabilization of the gemcitabine inside the host molecule, we carried out quantitative MESP, AIM analysis, non-covalent interaction analysis-reduced density gradient (NCI-RDG) analysis, and energy decomposition analysis (EDA) analysis.

Molecular electrostatic potential analysis
To corroborate the findings of ECT, we have computed the quantitative molecular electrostatic potential values for the gemcitabine, Q[n], and encapsulated complexes at (1) ECT = (ΔN max ) host − (ΔN max ) guest 0.001 electrons/Bohr 3 isodensity surface [55]. The MESP mapped isosurface for gemcitabine and Q[n] molecules is depicted in Fig. S5 and the inclusion complexes are shown Fig. 3. The electron-rich regions are represented in blue color and the deficient regions in red color and are assigned as V s,min and V s,max respectively.
In the guest molecule, we observed the V s,min on the oxygen atoms of carbonyl and hydroxyl groups. The highest V s,min was observed on the carbonyl group of the uracil unit. The V s,max was found on the hydrogen atoms of the amino group and hydroxyl groups. The 2′-OH has the highest value of 74.3 kcal mol −1 , while the 5′-OH has a value of 72.9. The amino hydrogen has 60.2 and 40.2 kcal mol −1 V s,min values. In the Q [6] host, the V s,min values are mainly localized on the carbonyl portal with a value of −102.7 kcal mol −1 . Besides, the inner cavity has a V s,min value of −51 kcal mol −1 . The V s,max were mainly localized on the methylene hydrogen atoms connecting the glycoluril units. In the Q [7] and Q [8] host, the V s,min on the carbonyl portal were −97.3 and 98.7 kcal mol −1 . Thus, with the increase in the size of the Q[n], only a minor variation is observed. Similarly, the V s,max also shows a small variation in size. However, the V s,min value observed on the inner cavity increases with the increase in Q[n] size.
It is anticipated that during the complex formation process, the noticed positive region in a molecule can interact with any negative region of the other molecule, giving rise to a directional interaction [56]. Hence, there could be an increase or decrease in the positive or negative values for the studied complexes. In the stable gemcitabine@Q [6] complex, the V s,max observed on the hydrogen atom of amino of gemcitabine turns to negative V s,max from 60.2 to −24.0 kcal mol −1 , while there is a substantial decrease in the V s,min value on both the carbonyl portal ends. Furthermore, the carbonyl portal on the end attached to the sugar ring of the host underwent a complete electron density shift, which is evident from the complete disappearance of blue regions on the carbonyl part of Q [6] host, while the V s,max observed on the methylene bridges shows decrease or increase in value depending on the location of the host molecule. In the case of gemcitabine@Q [7], although a similar trend has been observed, the charge transfer is more pronounced which is evident from the partial disappearance of red V s,max regions observed at the middle of the Q [7] unit. In the case of gemcitabine@Q [8] complex, a substantial change in the V s,max and V s,min values is observed on the cavity of the Q [8] unit; however, no appreciable change in the V s,min values of the carbonyl unit is observed. Thus, one can conclude that the MESP analysis shows electron delocalization is more pronounced in the gemcitabine@Q [7] complex, while in the gemcitabine@Q [8] complex, the complexation occurs at the inner cavity, which leads to a decrease in the electron delocalization compared to the other host-guest complexes.

Atom in molecule analysis
Atom in molecule analysis is a technique used to quantify the nature of the interaction between the molecules in a complex. The topological parameters such as total electron density ρ(r) and its Laplacian (∇ 2 ρ(r)) at the bond critical points provide the strength of bonds in the encapsulated complexes [57]. The computed total AIM parameters for the complexes are provided in Table 3, while their individual parameters are listed in Table S1-S3. The tomographs obtained for all the stable inclusion complexes are depicted in Fig. 4, in which the green dots represent the (3, −1) bond critical points, red dots (3, −1) ring critical points, and blue dots (+ 3, + 3) the cage critical points. The total number of intermolecular BCPs and their electron density were found to vary with the host molecule. The highest number of BCPs, ring critical points (RCPs), and cage critical points (CCPs) are found for the gemcitabine@Q [6] and least for the gemcitabine@Q [8] complex. The presence of the highest number of RCPs and CCPs in the gemcitabine@Q [6] complex shows the existence of steric crowding in the system which is a destabilizing factor for the complex.
The computed total electron density ρ(r) at the intermolecular BCPs have values in the range of 0.0005-0.0305 a.u. which are classified as weak bonding. The ρ(r) values of the hydrogen bonds were in the range of 0.0070-0.0257 a.u which is within the range proposed by Rozas et al. [58,59]. The highest ρ(r) values are found for the HOST C = O···H-N GUEST followed by HOST C = O···H-O GUEST . Besides the above interaction, we also noticed the presence of HOST O = C···F-C GUEST which is due to the halogen bonding [60]. Furthermore, we also noticed the presence of several N···N, O···O interactions in these complexes which are considered extremely weak. In addition, we observe the presence of one, two, and three new intramolecular bonds Fig. 3 Electrostatic potential maps of encapsulated complexes of a gemcitabine@Q [6], b gemcitabine@Q [7], and c gemcitabine@Q [8] superimposed on the isodensity surface of the structures (isovalue 0.001), computed at M05-2X-D3/6-311G(d,p) level Table 3 The average of QTAIM parameters (in a.u.) of intermolecular bonding between the gemcitabine drug molecule and cucurbituril Q[n] (n = 6, 7, and 8) hosts in gemcitabine molecules in the encapsulated complexes with host Q [6], Q [7], and Q [8] host molecules. In the gemcitabine@Q [8] complexes, we also noticed the presence of new intramolecular bonds between the carbonyl oxygen atoms indicating a high strain in the host Q [8] molecule. The ellipticity (ε) value is a measure of charge distribution, and thus the bond stability or weakness of the bond [61]. The ε value is > 1 for the gemcitabine@Q [6] complex which also supports its least stability. The sign of ∇ 2 (r) value is a useful measure to determine the nature of bonding as covalent or electrostatic and or closed-shell interaction [62]. The positive ∇ 2 (r) values observe at the intermolecular BCPs for all the complexes designate the interactions are electrostatic or closed-shell interactions. Besides, a positive value of H(ρ) designates the presence of non-covalent interaction at the BCP. Also, if the ratio of −G(ρ)/V(ρ) > 1, then the interaction can also support the existence of noncovalent interaction between the guest and host molecules. Thus, these complexes are stabilized by various non-covalent interactions including HBs and C···F interactions.

NCI-RDG analysis
NCI-RDG analysis helps to portray the various non-covalent interactions present in the complexes [63,64]. As these complexes are stabilized mainly by non-covalent interactions, NCI-RDG analysis is carried out on the most stable optimized structures. The 2D NCI analysis along with RDG isosurfaces for the encapsulated complexes is shown in Fig. 5.
In the 2D RDG plots, the strong hydrogen bonding interactions are shown in blue color (ρ > 0.01 a.u), while the green and red color spikes represent the van der Waals interaction and steric repulsion respectively. The sign(λ 2 )ρ values indicate that these encapsulated complexes have strong HBs, weak van der Waals, and steric interactions. From the isosurface NCI diagram, predominant steric interactions could be found in the gemcitabine@Q [6] complex which is evident from the red patches between the host and guest molecule. The NCI isosurface for gemcitabine encapsulated complexes with Q [7] and Q [8] host displays that the green patches are uniformly distributed in all directions and less amount of red patches providing them more stability compared to the gemcitabine@Q [6] complex. Thus, gemcitabine@Q [6] is destabilized by steric crowding, while van der Waals and hydrogen bonding are the main stabilizing factors in the formation of these complexes.

Energy decomposition analysis
To quantify the strength of intermolecular interactions between the gemcitabine and Q[n] host molecules in the encapsulated complexes, we used the EDA [65]. In EDA Fig. 4 Molecular topography analysis for the encapsulated complexes of a gemcitabine@Q [6], b gemcitabine@Q [7], and c gemcitabine@Q [8] Fig. 5 The gradient isosurfaces and the scatter diagrams for a gemcitabine@Q [6], b gemcitabine@Q [7], and c gemcitabine@Q [8] analysis, the total interaction energy is decomposed into Pauli's repulsion and attraction terms such as electrostatic, orbital, and dispersion energy using the Ziegler-Rauk Eq. (2) where ΔE Pauli is Pauling's repulsion energy and ΔV electrostatic , ΔE orbital , and ΔE disp are electrostatic interaction, orbital, and dispersion energies. ΔE Pauli is always a positive term owing to the destabilizing interactions linked with suitable antisymmetrization and renormalization of the Hartree product of the fragment wavefunctions and is accountable for steric repulsion. ΔE orbital is stabilizing component originating from the interactions of the occupied MOs of one fragment and the unoccupied MOs of the other fragment. ΔV electrostatic parallels the conventional electrostatic interaction between unperturbed charge distributions between the fragments. The calculated EDA values along with the percentage of contribution for the encapsulated complexes are listed in supporting information Table S4. A pictorial view of the same is provided in Fig. 6. The computed EDA results are quantitative and are not qualitative, as the dispersion interaction term needs explicit consideration of bulk water molecules, which is not present in the DFT and MD calculations [40]. Besides, recent studies have shown that cavitation energies can outperform dispersion interactions [47].
(2) ΔE int = ΔE pauli + ΔV electrostatic + ΔE orbital + ΔE disp It is apparent that the gemcitabine@Q [6] complex has the maximum Pauling's repulsive energy and is the least for the gemcitabine@Q [7]. In attractive terms, ΔV electrostatic is higher for the gemcitabine@Q [6], while for the other complexes, they are nearly identical. The orbital interactions are higher for the gemcitabine@Q [7] and least for the gemcitabineQ [6]. The dispersion energy is least for the gemcitabine@Q [6] and is largest for the gemcitabineQ [7] implying that dispersion plays a major role in stabilizing the complexes. Also, it can be seen that dispersion and electrostatic interactions in the gemcitabine@Q [7] complex have a nearly identical contribution. These results from the EDA analysis support the results obtained in the MESP analysis and NCI-RDG analysis. Thus, from the EDA outcomes, we can conclude that non-covalent and electrostatic interactions with partial covalent character alleviate the encapsulated complexes.

Conclusion
In this work, the complex formation capability of gemcitabine drug with host cucurbit[n]uril Q[n] (n = 6, 7, and 8) molecules in the gas and solution phase was carried out using density functional theory (DFT). The DFT study demonstrates Fig. 6 Polar graph of EDA results with contribution from Pauli's repulsion (black color) electrostatic attraction, orbital interaction, and van der Waals interaction. The attractive terms are provided in percentage, while the repulsion term is provided in kcal mol −1 that the most stable configuration is a fully encapsulated complex. In the gemcitabine@Q [6] and gemcitabine@Q [7] encapsulated systems, the gemcitabine amino -NH2 and the alcoholic group in the carbohydrate are bonding with the carbonyl units of Q[n]. In the case of gemcitabine@Q [8], the guest orientation differs and was parallel to the carbonyl portal. The computed deformation energy shows that the gemcitabine@Q [7] complex has the least deformation energy both on the host and guest molecules. The addition of sodium ions leads to the partial exclusion of the guest molecule. Furthermore, the formed complex loses the HBs between the gemcitabine and Q [7] and the sodium atoms lie close to the carbonyl portal of Q [7].
The computed binding energy of the complex was highest for the gemcitabine@Q [7] complex and least for the gemcitabine@Q [6] complex. Thermodynamic parameters computed for the complexation process exhibit high negative ΔG° and ΔH° values implying that the formation of the encapsulated complexes is spontaneous and is an enthalpydriven process. The entropy contribution for all the complexes was largely negative, implying that for the complexation process enthalpy plays a major role over entropy. Frontier molecular orbitals, namely HOMO and LUMO are located mainly on the gemcitabine uracil ring, before and after encapsulation formation, indicating that the encapsulation does not alter the electronic nature of the drug and the process occurs by pure physical adsorption. GRD indicate the charge transfer occurs from the guest to the host molecule and the charge transfer from the guest is nearly identical for all the host molecules. Thus, in addition to the charge transfer, other intermolecular interactions also play a vital role in stabilizing the guest inside the host molecule.
In the quantitative molecular electrostatic potentials of stable gemcitabine@Q [6] complex, the V s,max observed on the hydrogen atom of the amino group of gemcitabine turns to negative V s,max while there is a substantial decrease in the V s,min value on both the carbonyl portal ends. Furthermore, the carbonyl portal on the end attached to the sugar ring of the host underwent a complete electron density shift, which is evident from the complete disappearance of blue regions on the carbonyl part of Q [6] host. In the case of gemcitabine@Q [7], a more pronounced electron density shift was observed which is evident from the partial disappearance of red V s,max regions observed at the middle of the Q [7] unit. AIM topological analysis illustrates that these complexes are stabilized by various non-covalent interactions including HBs and C···F interactions. Besides the charge transfer, the computed topological parameters at the BCPs demonstrate the presence of non-covalent interactions between the host and guest molecules. The 2D RDG plots show the presence of strong HBs and weak van der Waals interaction and the presence of steric repulsion. The isosurface NCI diagram envisages predominant steric interactions in the gemcitabine@Q [6] complex which is evident from the red patches between the host and guest molecule. The NCI isosurface for gemcitabine encapsulated complexes with Q [7] and Q [8] host shows that the green patches are evenly distributed in all directions. EDA results demonstrate Pauling's repulsive energy which is a destabilizing factor to be predominant in the gemcitabine@Q [6] complex, while the orbital and dispersion energies stabilize the gemcitabine@Q [7] complex. Thus, from the above results, we can conclude that non-covalent and electrostatic interaction with partial covalent character stabilizes the gemcitabine encapsulated Q[n] complexes.