Heterofullerene MC59 (M = B, Si, Al) as Potential Carriers for Hydroxyurea Drug Delivery

As a representative nanomaterial, C60 and its derivatives have drawn much attention in the field of drug delivery over the past years, due to their unique geometric and electronic structures. Herein, the interactions of hydroxyurea (HU) drug with the pristine C60 and heterofullerene MC59 (M = B, Si, Al) were investigated using the density functional theory calculations. The geometric and electronic properties in terms of adsorption configuration, adsorption energy, Hirshfeld charge, frontier molecular orbitals, and charge density difference are calculated. In contrast to pristine C60, it is found that HU molecule is chemisorbed on the BC59, SiC59, and AlC59 molecules with moderate adsorption energy and apparent charge transfer. Therefore, heterofullerene BC59, SiC59, and AlC59 are expected to be promising carriers for hydroxyurea drug delivery.


Introduction
Nanomaterials, such as C 60 and its derivatives, have become increasingly important in gas sensors and biomedicine, especially in the field of drug delivery [1][2][3][4][5][6][7][8]. As the routine methods of drug administration, which suffered from the problems of lacking target specificity and unavoidable side effects, are in great need of improvement [9], much efforts have been done to develop new drug delivery systems based on nanomaterials. Nowadays, drug carriers based on nanomaterials can be used to control the release of a drug into systemic circulation, which could be achieved by triggered release at the target site by stimulus, such as changes in pH, application of heat, and activation by light [9].
In recent years, fullerenes have been regarded as one of the most exciting materials for developing new drug delivery systems, due to their unique properties including small dimensions, high surface-to-volume ratios, ease of functionalization, and biocompatibility [10][11][12]. For example, Wang et al. have successfully designed a nanoplatform based on a nucleus-like fullerene core and realized the precisely controlled delivery of theranostic agents both in vitro and in vivo [13]. Shi et al. found that PEI-derivatized fullerene could be used as drug delivery vehicles to achieve suppression of tumor growth without toxic effects on normal tissues [14].
Hydroxyurea (HU), also known as Hydroxycarbamide, is an effective anti-cancer drug widely used in treating head and neck cancer, breast cancer, chronic myelocytic leukemia, and so on [33,34]. However, studies have shown that HU drug can cause some side effects, such as rashes, leukopenia, and leg ulcers during treatment [35,36], which significantly limits the application of HU in the treatment of cancer. Therefore, in order to reduce side effects and enhance the effectiveness, it is necessary to find a suitable drug delivery vehicle for HU drug [37,38]. Considering that the experimental procedure is expensive and time-consuming, we performed theoretical calculations with the aim to find potential carriers for HU drug delivery, which could offer more options to the experimental work.
In this work, we employ density functional theory calculations to accurately describe the adsorption of HU drug on the pristine C 60 and heterofullerene MC 59 (M = B, Si, Al), in order to reveal some clues for a drug delivery vehicle. In the following section, we present the details of the computational methods, which are then followed by the results and discussion. Finally, the general overview and conclusions are given.

Computational Details
In this work, all the geometry optimization and energy calculations were performed using the spin-polarized density functional theory (DFT) implemented in the DMol 3 program package [39,40]. The gradient-corrected (GGA) exchange-correlation functional formulated by Perdue, Burke, and Ernzerhof (PBE) [41] along with the double numerical basis sets including d-polarization functions (i.e., the DNP set) are chosen here. In order to better describe the weak interaction between fullerene and drug molecules, we introduce a long-range dispersion interaction correction term produced by Grimme during the calculation [42]. In the self-consistent field calculations, the convergence criterion was set to 10 −6 a.u. for energy and electron density. The geometries are fully optimized with no restrictions. In geometrical optimizations, we set the convergence criterion of 10 −3 a.u. on the gradient and displacement and 10 −5 a.u. on the total energy. To make sure the resultant structures belong to the real local minima, the normal-mode vibrational analysis was applied. No imaginary frequency was observed through vibration frequency analysis. The charge transfer from the HU molecule to fullerene is analyzed based on Hirshfeld charge analysis [43], which is based directly on the electron density as a function of space.
The adsorption energy (E ads ) was defined as follows: where E complex , E HU , and E fullerene or heterofullerene are the total energy of the HU adsorbed on C 60 or B-, Si-, and Al-C 59 fullerenes, the energy of the independent HU drug, and the energy of the C 60 and B-, Si-, and Al-C 59 fullerenes, respectively. In order to understand the degree of charge redistribution between the drug molecule and fullerene in more detail, charge density difference (CDD, ∆) maps were calculated and defined as follows: where ρ complex , ρ HU , and ρ fullerene or heterofullerene are the electron density of the HU adsorbed on C 60 or B-, Si-, and Al-C 59 fullerenes, the electron density of the independent HU drug, and the electron density of C 60 and B-, Si-, and Al-C 59 fullerenes, respectively. In the CDD calculation, we put the complex into a periodic 25 × 25 × 25 Å cubic box, and the k-sampling was restricted to the Γ point. The other calculation parameters are the same as described above.

The Hydroxyurea Molecule Characterizations
The optimized geometry, molecular electrostatic potential (MEP) plot, highest-occupied molecular orbital (HOMO) and lowest-unoccupied molecular orbital (LUMO) plots of hydroxyurea (HU) are presented in Figure 1. The MEP plot shows that the O(2) atom of HU has the highest negative electrostatic potential and, thus, could be considered as the active site when interacting with fullerene molecule [19][20][21]. Although the two N atoms also have slightly negative electrostatic potential, N atom sites probably cannot be the adsorption sites due to the steric effect [22]. Moreover, the highest density regions of HOMO and LUMO profiles are mainly localized on the O(2) atom and the H atom adjacent to O(1) atom, respectively, which is consistent with the blue and red regions on the MEP plot.

The Adsorption of HU on the Pristine C 60
C 60 is composed of 20 hexagonal and 12 pentagonal rings, and C-C bonds can be classified into two types. The first one is shared between two hexagons and referred as  bond, and the second one is shared by one hexagon and one pentagon rings and named as  bond. As shown in Figure 2, the calculated bond lengths in C 60 are 1.451 and 1.399 Å, essentially in good agreement with experimental values [44]. There are several possible adsorption configurations when HU molecule is attached to C 60 . All the high symmetry adsorption sites, including the five-membered ring, the six-membered ring, the C-C bond and C atom top, have been considered as the active sites of the C 60 . As shown in Figure 2, the most stable configuration is the H atom of HU adsorption on the top of the C atom with a distance of 2.24 Å. Upon adsorption, both the bond lengths around the adsorption C atom of C 60 and HU molecule are nearly unchanged. In addition, the adsorption energy (Table 1) of the pristine C 60 and HU systems is −0.271 eV. The relatively low adsorption energy and the almost unaltered geometric parameters of pristine C 60 and the HU indicate that the interaction between them mainly occurs via noncovalent bond [45].  In order to better understand the interaction between pristine C 60 and HU drug molecule, we further analyzed the charge transfer based on Hirshfeld charge analysis. The charge transfer from HU to pristine C 60 is −0.047 e (Table 1), which is relatively small. Moreover, as can be seen from Table 1, the HOMO and LUMO energy levels of HU-C 60 complex have only slightly changed compared to the initial C60, resulting in the HOMO-LUMO energy gap (Eg) also reduced very slightly, which suggests the adsorption of HU molecules has little effect on the electronic properties of C 60 [19][20][21][22]. Figure 2c shows the charge density difference (CDD) plot of HU-C 60 with the isosurface value of 0.03 a.u. It is apparent from this figure that few charge density difference in the region between HU and C 60 can be observed. Besides the lower adsorption energy and charge transfer, and slightly changed energy gap, the insignificant charge density difference further confirms the weak interaction nature between HU and C 60 [17,18]. Therefore, pristine C 60 is not a proper drug carrier for HU. To enhance the reactivity of C 60 , we explored the effect of replacing a C atom of C 60 with B, Si, and Al atom on the adsorption properties of HU drug.

The Adsorption of HU on the Doped C 60
The optimized structures of the doped C 60 by B, Si, and Al atom are presented in Figure 3. Due to the larger covalent radius, all the dopant atoms caused certain deformations at the point where they were inserted in C 60 cage. Compared with pristine fullerenes, the geometry of BC 59 undergoes a slight deformation, and the corresponding (6-5) and (6-6) B-C bond lengths are stretched to 1.545 and 1.492 Å, respectively. In SiC 59 , the formed (6-5) and (6-6) Si-C bond lengths are 1.843 and 1.793 Å. Compared with BC 59 and SiC 59 , the structural deformation caused by Al doping becomes more pronounced, since the Al atom is larger in size. Meanwhile, the M-doping also significantly modifies the electronic structures of the C 60 . It can be seen from Table 1 that the HOMO-LUMO energy gap (E g ) of the doped system is significantly decreased, especially for the BC 59 and AlC 59 . In the framework of conceptual density functional theory [46], the lower E g means higher reactivity. Therefore, doping by B, Si and Al atoms has improved the reactivity of C 60 remarkably. In SiC 59 , the decrease of E g is less significant than the case of BC 59 and AlC 59 , which could be explained by the Si atom is valence isoelectronic with C atom [7].
According to the above analysis, the doping with B, Si, and Al atoms properly modify the structure and electronic properties of C 60 , which may make it possible for the delivery of HU molecule. The computed MEP plots of BC 59 , SiC 59 , and AlC 59 are shown in Figure 3. The MEP plot, which is based on the distribution of charge density, is an efficient approach for the visualization of the reactive sites on the molecule. It is apparent from Figure 3 that the B, Si, and Al sites have the highest positive electrostatic potential in BC 59 , SiC 59 , and AlC 59 , respectively. So, the B, Si, and Al atoms act as the electrophilic sites. As mentioned above, the O(2) atom in HU molecule act as the nucleophilic sites. Therefore, we can predict that the HU molecule should absorb from its O(2) atom to the doping site of the MC 59 .
The optimized geometries of the HU-MC 59 complexes are shown in Figure 4. The most stable structures are HU molecule interact with the B, Si, and Al atoms (electrophilic sites) of MC 59 through its O(2) atom (nucleophilic site), which is in agreement with the MEP prediction. It can be seen from Table 1 that the adsorption energies of the HU-BC 59 , HU-SiC 59 , and HU-AlC 59 complexes are −1.122, −1.567, and −2.174 eV, respectively. A negative value of the adsorption energy indicates that the adsorption process is exothermic and geometrically stable. Generally, the energetic threshold separating the adsorption energy of physisorption and chemisorption is about 0.5 eV per adsorbed species [47]. Due to the relatively more considerable adsorption energy, the interaction between the HU and the MC 59 should be chemisorption [5][6][7][8]. It can be seen that the HU-AlC 59 mixture has the largest adsorption energy and is the most stable one among the three complexes.  Figure 2) for all. In the CDD plots, the blue and yellow colors correspond to the electron density gain and loss regions, respectively.
To further understand the interaction between HU drug molecule and MC 59 , we calculated the charge transfer based on Hirshfeld charge analysis. From Table 1, it can be seen that the amount of charge transfer of HU-MC 59 complex (0.373, 0.429, and 0.346 e) is much higher than HU-C 60 (−0.047 e). The CDD plots of HU-MC 59 are also calculated and shown in Figure 4. Different to the case of HU-C 60 (Figure 2), here we can observe obvious charge differences between the HU and MC 59 .
Furthermore, as shown in Figures 5 and 6, we have calculated the frontier molecular orbitals of the MC 59 molecules and HU-MC 59 complexes. It can be seen from Figure 5 that in MC 59 the HOMO and LUMO are mainly distributed around the M-dopant with different orbital shapes, which is in agreement with previous reports [21,22]. However, after the HU molecule adsorbed on the MC 59 (Figure 6), the original charge density above the M-dopant atom disappeared for both HOMO and LUMO, which indicates the relatively stronger interaction between HU and MC 59 [7].  Overall, the moderate adsorption energy, moderate charge transfer and the obvious charge density difference confirm the chemisorption nature of HU on MC 59 . These findings suggest that heterofullerene BC 59 , SiC 59 , and AlC 59 molecules can be used as the potential drug carrier for HU drug.
From the perspective of theoretical calculations, research strategies on the potential applications of heterofullerene in drug delivery [19,20,23] and drug detection [16,21,22] are similar. Hence, we further analyzed the possibility of MC 59 for HU drug detection.
We utilize the relation between electronic conductivity and band gap of an intrinsic semiconductor [48]: where σ is the electric conductivity, E g is the band gap, k is known as the Boltzmann's constant, and T is the thermodynamic temperature. It can be seen from Formula 3 that the electric conductivity is exponentially related to E g at a given temperature. Thus, a smaller change in E g results a larger variation in the electric conductivity. Although this formula is derived for intrinsic semiconductor, it has also been widely employed to study the potential applications of nanomaterials in gas sensing [5][6][7][8]49,50] and drug detection [16,22]. Take the Al-doped ZnO nanostructure (AZO) for CO chemical sensors for example [51,52]. The experimental results show that the electronic conductivity of the AZO film increases significantly upon CO adsorption [51], which can be understood well with the reduction of the HOMO-LUMO gap of AZO molecule after the adsorption of CO [52]. The validity of Formula 3 for nanomaterials is mainly due to the relatively weak interaction between nanomaterial molecules, which makes the nanomaterial films mimic the electronic properties of the single molecule to some extent. After the HU adsorption, the HOMO-LUMO energy gap (E g ) of HU-MC 59 complexes change obviously (in the range of 13.4% to 33.7% Table 1), which would result a larger variation in the electric conductivity of MC 59 film. Therefore, the HU molecule could be detected by MC 59 (M = B, Si, Al) theoretically.
At the end of this article, we briefly discuss the toxicity of fullerene and its derivatives. Fullerene and its derivatives, which have been extensively studied in biomedicine, have shown excellent potential applications in many fields, such as cancer therapy and photodynamic therapy [10][11][12][13][14]. However, as stated in some review articles, fullerene and its derivatives' presence or absence of toxicity remains a controversial issue [53,54]. Some research results show that fullerene derivatives are toxic. On the other hand, a series of tests show that fullerene derivatives are non-toxic, especially at low concentrations [53,55]. For example, the recent experimental work by Liu et al. [56] has shown that fullerene derivative Gd@C 82 (OH) 22 can be used as a non-toxic breast cancer stem cell-specific inhibitor. Therefore, it should be pointed out that our theoretical calculation results here offer a potential carrier for HU drug, and further investigations on toxicity and clinical trials are indispensable before application.

Conclusions
In summary, with the aim to find a promising nanocarrier for HU drug delivery, we have investigated the interactions between HU drug and pristine C 60 and heterofullerene MC 59 (M = B, Si, Al) using density functional theory calculations in terms of adsorption geometries, adsorption energies, charge transfer, and electronic properties. It is found that the HU molecule is physisorbed on the pristine C 60 with lower adsorption energy and negligible charge transfer. So, pristine C 60 cannot act as the carrier for HU drug delivery. On the contrary, the HU molecule is chemisorbed on the BC 59 , SiC 59 , and AlC 59 molecules with strong binding and can lead to finite charge transfer. Thus, it is suggested that heterofullerene BC 59 , SiC 59 , and AlC 59 are potential candidates for HU drug delivery.
Author Contributions: P.W. and J.Z. proposed the project and analyzed the simulation results. G.Y. and P.W. contributed to the DFT simulations. X.Z., Y.D., D.C., and J.Z. performed the analysis of the data and provided some revision of the manuscript. All authors have read and agreed to the published version of the manuscript.