Surface modification of graphene and fullerene with Sulfur (S), Selenium (Se), and Oxygen (O): DFT Simulation for enhanced zidovudine delivery in HIV treatment

Abstract

HIV is one of the most threatening health conditions with a highly increasing rate, affecting millions of people globally, and from its time of discovery until now, its potential cure cannot be explicitly defined. This challenge of having no/low effective drugs for the subjected virus has called for serious attention in the scientific world of virus disease therapeutics. Most of these drugs yields low effectiveness due to poor delivery; hence, there is a need for novel engineering methods for efficient delivery. In this study, two nanomaterilas (graphene; GP, and fullerene; C60) were modelled and investigated with sulfur (S), selenium (Se), and oxygen (O) atoms, to facilitate the delivery of zidovudine (ZVD). This investigation was computationally investigated using the density functional theory (DFT), calculated at B3LYP functional and Gd3bj/Def2svp level of theory. Results from the frontier molecular orbital (FMO), revealed that the GP/C60_S_ZVD complex calculated the least energy gap of 0.668 eV, thus suggesting a favourable interactions. The study of adsorption energy revealed chemisorption among all the interacting complexes wherein GP/C60_S_ZVD complex (-1.59949 eV) was highlighted as the most interacting system, thereby proving its potential for the delivery of ZVD. The outcome of this research urges that a combination of GP and C60 modified with chalcogen particularly, O, S, and Se can aid in facilitating the delivery of zidovudine.

Introduction

Human immunodeficiency virus/acquired immunodeficiency syndrome (HIV/AIDS) presents a formidable global health challenge with far-reaching implications. This virus targets the immune system, undermining the body’s ability to combat infections and diseases effectively. Consequently, individuals infected with HIV are highly susceptible to infections and diseases [1, 2]. It is indisputable that HIV is a significant health concern that affects millions of people worldwide and a definitive cure for HIV has not yet been identified [3,4,5,6,7,8]. The challenge has attracted significant attention in the scientific community of drug manufacturing [4]. Antiretroviral drugs, such as zidovudine (ZDV), lamivudine, stavudine, and dideoxycytidine, have emerged as effective agents for managing HIV infections [9,10,11,12]. However, challenges persist in their effective delivery to the target site within the host immune system. Issues such as limited penetration, adverse side effects, and metabolic concerns underscore the need for innovative drug delivery systems [13]. However, zidovudine, also known as azidothymidine (AZT), is a nucleoside analog reverse transcriptase inhibitor (NRTI) that plays a crucial role in the management of human immunodeficiency virus (HIV) infection. It was the first antiretroviral drug approved by the FDA in 1987 and remains a cornerstone of combination therapy for HIV-1 infection [14]. The mechanism of action of zidovudine involves the inhibition of reverse transcription, a critical step in the HIV replication cycle. Zidovudine is phosphorylated to its active metabolite, zidovudine triphosphate, which competes with thymidine triphosphate for incorporation into viral DNA. This leads to the termination of DNA synthesis, thereby inhibiting viral replication [15]. Despite its efficacy, zidovudine delivery to target HIV-infected T-cells remains a challenge. The drug has poor bioavailability, rapid clearance, and limited penetration into lymphoid tissues, where HIV reservoirs exist [16].

To improve zidovudine delivery, recent advancements in nanotechnology and computational chemistry approaches offer promising solutions to address these challenges [17,18,19]. Many studies have explored various strategies, including developing nano formulations and drug delivery systems. For instance, Petlin et al. (2017) developed a zidovudine-loaded nanoparticle formulation that demonstrated enhanced bioavailability and antiretroviral activity compared to the free drug [20]. Another study by Ilinskaya et al. (2019) designed a zidovudine-conjugated nanoparticle that targeted HIV-infected T-cells, resulting in improved drug delivery and anti-HIV activity [21]. Amin and Boateng (2020) discussed the use of ligand-targeted liposomes and nanoparticles for delivering zidovudine to HIV-infected cells [22]. Another review by Pradhan et al. (2021) on Drug Delivery highlighted the potential of nanotechnology-based approaches in enhancing zidovudine delivery and reducing toxicity [23]. Also, Dahmane et al. (2013), Joshy et al. (2018), and Gupta et al. (2019) demonstrated the feasibility and efficacy of utilizing various nanoparticles, including chitosan, PF-68-coated alginate conjugates, and fullerene-like nanocages for drug delivery [24,25,26]. Their outcome unveils the essentials of nanomaterials in drug delivery applications. Therefore, in this research, we focused on improving the delivery of zidovudine by engineering graphene sheet and fullerene C60 together and also doped the fullerene structure with chalcogens such as oxygen, selenium, and sulfur. Integrating chalcogen elements (oxygen, sulfur, and selenium) into the fullerene structure can enhance the material’s biocompatibility and drug-loading capacity [27,28,29]. The chalcogen-doped nanocomposite can exhibit a high surface area and excellent conductivity that can be a suitable platform for drug carriers and controlled release [30,31,32]. This composite’s unique physicochemical properties contribute to improved drug stability, prolonged circulation time, and enhanced cellular uptake, ultimately optimizing the therapeutic efficacy of a drug. Using fullerene C60 as a nanovesicle can be attributed to the unique biocompatibility essential for drug delivery, as reported by Guan et al. and Lin et al. (2016 & 2017) [33, 34]. The unique electronic and chemical properties of oxygen, selenium, and sulfur have garnered attention [35, 36], and these properties make them ideal candidates for incorporation into fullerene structures aimed at improving the surface functionality and opening paths for tailoring their interaction with biological entities as well as enabling specific recognition and uptake by HIV-infected cells [37, 38]. The potential to deliver zidovudine (ZDV) directly to the affected cell holds promise for maximum efficacy while minimizing off-target effects and side effects.

Our study employed an advanced computational density functional theory approach encompassing electronic properties studies (frontier molecular orbital, natural bond orbital, and density of state) to elucidate the chemical reactivity and kinetic stability mechanisms between nanocomposites and ZDV. The study also utilized topological analysis to effectively assess natural interactions and interatomic interactions. The noncovalent interaction (NCI) and quantum of atom in molecule (QTAIM) approaches were instrumental in obtaining these results. To determine whether the drugs will be chemically or physically adsorbed, an adsorption study and other objectives for desorption mechanism were achieved successfully via density functional theory (DFT) calculations. Continued research in this field is essential to refine and optimize these nanocomposites for enhanced antiretroviral drug delivery. Ultimately, this study revealed promising outcomes for improving treatment plans for individuals with HIV/AIDS.

Method

DFT approach

This computation was carried out employing first-principles density functional theory (DFT) simulations, calculated at B3LYP/D3(BJ)/Def2svp levels of theory. The optimization of the combined graphene sheet and fullerene C60 nanomaterials doped with chalcogen groups was performed utilizing B3LYP functional and the D3 (BJ) dispersion correction in conjunction with the Def2svp basis set, a choice given the need for a balanced compromise between accuracy and computational efficiency in describing the intricate interactions within these complex nanocomposites. This method combination has been extensively validated in previous studies [39, 40] and is particularly well-suited for capturing the subtleties of non-covalent interactions and chemical bonding in large-scale systems, thereby providing a robust framework for elucidating the structural and electronic properties of these materials. The structural optimization and electronic properties of GP/C60 were analysed to understand the nature of the coordination between selenium, oxygen, and sulfur dopants and the graphitic carbon matrix, which was carried out using Gaussian 16 [41] and Gaussian view 6.6.0.16 [42]. This approach was employed to tailor accurate theoretical calculations. For proper investigation of the electronic, structural, and topological properties for effective delivery of zidovudine, an allusive objective was carried out. The reactivity and molecular stability of the surface doped with ZVD were considered vital for understanding the delivery dynamics; thus, frontier molecular orbital (FMO) analysis has been a powerful computational tool employed to investigate various chemical processes, including the possibility of electron transfer from individual doped and co-doped surfaces to drugs, by considering the energy gap obtained from the value of the HOMO and LUMO. The iso-surface of the frontier molecular orbital was visualized using Chemcraft software version 1.6 [43]. Similarly, other investigations were carried out using Multiwfn package 3.7 [44] for the quantum theory of atom-in-molecule (QTAIM) studies, which enabled the study of the intermolecular interactions occurring in the system, and the electronic distribution study was performed using the DOS approach also utilizing the Multiwfn package. Furthermore, the study of the non-covalent interaction (NCI) was carried out using the visual dynamic simulation (VMD) software package 1.9.4 [45]. Furthermore, the interactions between zidovudine and GP@C60_Se, GP@C60_S, GP@C60_O, and GP@C60_O_S_Se were explored to elucidate the drug-loading mechanism and release dynamics. The impact of pH changes on drug release was pivotal in this study, and computational simulations were conducted to elucidate these effects. The simulations utilized density functional theory (DFT) calculations implemented in Gaussian 16 software to predict the protonation states and electrostatic interactions of zidovudine and the GP/C60 carrier at different pH levels. This approach allowed for the modelling of how pH variations influence the solubility, ionization, and binding affinities within the drug delivery system. Additionally, to gain an in-depth understanding of the chemistry of 2D nanostructures, structural studies and analyses of binding energy, vibrational stability, and UV‒vis excitation was carried out, providing perspectives on the design and optimization of GP/C60 nanostructures for efficient and controlled therapeutic delivery of zidovudine and advancing antiretroviral drug delivery systems.

Results and discussion

Chemistry of 2D nanostructures

Structural studies

This research comprehensively examined the structural properties of a nanocomplex and its doped counterpart (see Fig. 1) as well as their interaction with zidovudine (see Fig. 2) by employing density functional theory (DFT) to elucidate changes in bond lengths. This investigation focused on understanding the implications of drug interactions on the potency, likely interactions, and stability of the nano complexes, which will determine their suitability for efficient delivery [46]. The bond length, which is the distance between the nuclei of atoms within a chemical bond, is a crucial parameter for this analysis [47]. This study specifically explored the effects of the GP/C60_O_ZVD nanocomplex and revealed a slight increase in the bond length before drug interaction, with minimal variations post-absorption. The differences in bond length, which are methodically recorded in Table S1, revealed a sequential change in the number of atomic nuclei in the nanocomplexes. For instance, in the case of O96—C26, there is a marginal difference of 0.00119 Å, O96—C25 exhibits a variation of 0.00314 Å, and C25—C49 displays a difference of 0.0006 Å. The overall change in bond length across all atomic nuclei is approximately 0.001 Å, suggesting that the interactions do not significantly decrease the bond lengths. Furthermore, this study was extended to the analysis of other nanocomplexes, such as GP/C60_O_S_Se_ZVD, GP/C60_S_ZVD, and GP/C60_Se_ZVD. GP/C60_O_S_Se_ZVD exhibited diverse bond lengths, with Se96—C25 displaying a greater bond length than the other atomic pairs. Conversely, GP/C60_S_ZVD reveals a substantial variation in bond lengths before and after drug interactions, particularly with a 0.03 Å change in some bonds and a lower 0.009 Å variation in carbon‒carbon bonds. Also, GP/C60_Se_ZVD demonstrated the greatest bond length, specifically for selenium-carbon bonds. However, the observed trend of decreasing bond lengths, as indicated by GP/C60_Se_ZVD > GP/C60_S_ZVD > GP/C60_O_S_Se_ZVD > GP/C60_O_ZVD, suggests a correlation with enhanced drug delivery potency. Because the decrease in bond length as observed across the systems depicted a sensitivity behaviour and also implies that the surfaces become more conjugated which is vital characteristics of drug delivery systems. Therefore, these systems can be beneficial in the development of drug delivery system.

Fig. 1

Structural representation of the investigated systems before and after interaction with ZVD

Fig. 2

3D Illustration of the drug target (zidovudine)

Vibrational stability

Vibrational analysis was carried out to determine the composition of the complexes in terms of the number of functional groups present [48]. This analysis identifies specific functional groups within molecules based on their characteristic absorption bands and this allow us to infer the types of ligands and metal ions present in the complex. The shifts in vibrational frequencies upon complex formation provide direct evidence of interactions between the metal ion and ligands thereby confirming the formation of the complex and indicating the nature of metal-ligand bonding (e.g., coordination bonds) [49]. Also, the intensity of absorption bands in the IR spectrum correlates with bond strengths; stronger bonds, such as metal-ligand bonds, typically exhibit more intense absorption peaks. This information helps assess the stability of the complex under various conditions [50]. Therefore, vibrational analysis through IR spectroscopy is essential for characterizing complex composition and stability, providing crucial molecular-level information for advancing our understanding of coordination chemistry and materials science. The frequencies are recorded for both surfaces after their interaction with ZVD to determine the stability of the nanocomposites. The CH functional group is predominant and appears in both symmetric and asymmetric patterns, as shown in Table S2. Functional groups such as OH, N₃, NH, and CO are only present in the interactions due to the effect of zidovudine. For the CH symmetric functional group, the frequencies before interaction ranged between 3423.64 cm^− 1 and 3191.82 cm^− 1. All the compounds exhibited an exact peak at approximately 3191 cm^− 1 except for GP/C60_O_S_Se_ZVD, which differs greatly from the peak at 3423.64 cm^− 1. After the interaction, the symmetric CH functional group exhibited an increase in the range of 3192.60 cm^− 1 to 3598.01 cm^− 1. In this case, GP/C60_O_S_Se_ZVD decreases slightly, while the other compounds increase. The asymmetric CH functional group also falls in a similar range as the symmetric group, and GP/C60_O_S_Se_ZVD still shows the same deviation. The functional group CO is recorded for only GP/C60_O_ZVD and GP/C60_Se_ZVD because the other two compounds did not experience vibrations at the CO functional groups. For the surfaces, only GP/C60_O_ZVD experiences a CO vibration since oxygen is part of its components at both the surface and interaction levels. For GP/C60_O_ZVD, the frequency before the interaction was recorded to be within the range of 1600 cm-1. After the interaction, the CO functional frequency ranged from 1800 cm^− 1 for GP/C60_O_ZVD to 3175.77 cm^− 1 for GP/C60_Se_ZVD, which is much greater than the normal range for the CO functional group [49, 50], indicating high stability of the CO functional group of the GP/C60_Se_ZVD complexes. The N₃ vibrational frequency is only recorded in GP/C60_O_ZVD and GP/C60_Se_ZVD and falls in the range of 1700 cm^− 1 for both compounds.

UV‒vis excitation analysis

The major aim of this analysis was to determine the behaviour of the complexes when the absorption of ultraviolet rays occurs at a particular wavelength. In Table 1, parameters such as the stabilization energy, excitation type, percentage transitions between the highest occupied molecular orbital and the lowest unoccupied molecular orbital, and the strength of the oscillators in each excitation type are presented. It has already been established from the energy of photons that energy increases with increasing frequency and decreasing wavelength [51], which is represented by the following equation:

$$\:E=h\times\:\frac{v}{\lambda\:}$$

where h is Planck’s constant (h = 6.626 × 10^− 34 Js^− 1), v is the speed of light (v = 3.0 × 108 ms^− 1) and λ is the wavelength at which the system absorbs ultraviolet radiation [52, 53]. The interactions are discussed for three excitation types: S₀-S_1, S₀-S₂, and S₀-S₃. For this study, as presented in Table 1, the stabilization energy increased in the following order: GP/C₆₀_O_ZVD > GP/C₆₀_S_ZVD > GP/C₆₀_Se_ZVD > GP/C₆₀_O_S_Se_ZVD. However, a slight difference is observed in GP@C₆₀_O_S_Se_ZVD, where the first excitation type S₀-S₁ is much lower by 0.763 eV compared to S₀-S₁ of GP/C₆₀_O_ZVD, with an energy of 0.885 eV, which is second to the former. Comparing the S₀-S₃ excitation type for the interactions, GP/C₆₀_O_S_Se_ZVD records the highest energy of 1.358 eV on the table, followed by GP/C₆₀_Se_ZVD with an energy of 1.159 eV in the S₀-S₃ type. The wavelength generally decreases from S₀-S₁ to S₀-S₃. This is due to the increase in energy in the system, as seen in the formula above. The wavelengths are generally longer in S₀-S₁ and shorter in S₀-S₃. For GP/C₆₀_O_ZVD, the energy decreases between 1147.13 and 1401.7 nm, that of GP/C₆₀_S_ZVD decreases between 719.13 and 983.77 nm, that of GP/C₆₀_Se_ZVD decreases between 1069.59 and 1222.35 nm, and that of GP/C₆₀_O_S_Se_ZVD decreases between 912.97 and 1624.60 nm. Comparing these interactions, we see that GP/C₆₀_O_ZVD records the largest wavelengths and smallest energy stabilization energy. This implies that the system is less stable and more reactive since its light absorption rates are so intense judging from its long wavelengths. The percentage contribution of each excitation type is recorded. GP/C₆₀_O_ZVD at S₀-S₁ with a H + 1-> L-1 transition indicates that the absorption occurs slightly below the HOMO and slightly above the LUMO with a percentage of 149.6%, and at S₀-S₂ with a H-> L + 1 transition occurs exactly between the HOMO and slightly below the LUMO, which still indicates that this particular excitation occurs at the HOMO with a contribution of 115.1%. The last excitation type, S₀-S_3, occurs directly between the HOMO and the LUMO with a contribution of 140.4%. The second interaction, GP/C₆₀_S_ZVD at the first excitation type S₀-S_1, with a transition from the HOMO to the LUMO of 140.7%, the second excitation type with a transition between the HOMO and slightly below the LUMO, which is still the HOMO with a contribution of 140.6%, and the third excitation type with a transition between the HOMO and two orbitals below the LUMO, which is one orbital below the HOMO with a contribution of 139.3%. Similar transitions are recorded for GP/C₆₀_Se_ZVD and GP/C₆₀_O_S_Se_ZVD. The percentage contribution ranges from 96.1% for GP/C₆₀_O_S_Se_ZVD S₀-S₃ to 164.9% for GP/C₆₀_O_S_Se_ZVD S₀-S₁. The oscillator strength (f) is a dimensionless quantity that ranges from 0 to 1 and provides information on the probability of the transition of electrons. From the table, GP/C₆₀_O_S_Se_ZVD records the highest transition of electrons since it possesses the highest oscillator strength of approximately 1.07 × 10^− 2 in its second excitation state. The oscillation strength of GP/C₆₀_Se_ZVD was 0.0, which indicates that no electrons were transferred in the third excitation type. However, the oscillator strength (f) is related to the dipole moment (µ) by the equation below:

f ∝ µ².

This means that as the oscillator strength increases, the dipole moment also increases, indicating a greater separation of charge and thus increased polarity. In the context of this study, the increased oscillator strength and corresponding dipole moment enhance the electron transfer process. A larger dipole moment facilitates the transfer of electrons between molecules as it creates a greater electric field gradient. This in turn, increases the probability of electron transfer, which is crucial for understanding the UV-Vis absorption properties of the system. As the oscillator strength increases, the dipole moment increases, and the polarity of the system increases. It was observed that GP/C₆₀_O_S_Se_ZVD and GP/C₆₀_S_ZVD are more polar than GP/C₆₀_O_ZVD and GP/C₆₀_Se_ZVD, suggesting to be more conjugated.

Table 1 UV‒vis analysis results for all the studied systems showing the following parameters: excitation type, energy (eV), wavelength (nm), percentage contribution (%), oscillator strength (f) and transition

Electronic property investigation

Chemical reactivity descriptors

The difference between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), known as the energy gap, is a crucial indicator for understanding the chemical complex reactivity, electronic structure, and stability of the complex [54]. The chemical and physical properties of a complex are significantly influenced by its molecular orbital, which is a function that characterizes the frontier molecular orbital (FMO) theory, as reported in previous studies. In this study, the energy gap is shown in Table 2 and illustrated in Fig. 3. The complex shows a variation in the surfaces and interacting complex, where the bear surfaces possesses a lower energy gap compared to the complex with the medication. As shown in the table, GP/C60_O_S_Se_ZVD had an energy gap of 1.276 eV before interaction and 1.303 eV after interaction. The GP/C60_O_ZVD complex had an energy gap of 1.699 eV before interaction and 1.778 eV after interaction, and GP/C60_Se_ZVD exhibited an energy gap of 1.959 eV before interaction and 2.083 eV after interaction. The energies of the GP/C60_O_S_Se_ZVD, GP/C60_O_ZVD, and GP/C60_Se_ZVD complexes increase in the same way after the interaction. However, for the complex GP/C60_S_ZVD, which has an energy gap of 1.895 eV before and 0.668 eV after interaction, a different trend is observed, which also gives it a greater chance of delivering the drug effectively than the other complexes. Chemical hardness and softness play a crucial role in confirming the results, which explains that the ability of a chemical substance to transfer electron charge within its immediate molecule depends on its chemical hardness and softness [55, 56]. Chemical hardness also influences the stability and reactivity of a complex, depending on the hardness of the complex [56]. As shown in the table, GP/C60_O_ZVD exhibited a high chemical hardness of 0.889 eV, and GP/C60_S_ZVD had the lowest chemical hardness of 0.334 eV. To buttress the finding, the complex GP/C60_S_ZVD is observed to have a very high chemical softness of 1.495 eV, suggesting the relative instability and reactivity of the complex. Upon interaction/complexing with the drug separately, a decrease in softness is observed, following this trend for the complexes. GP/C60_O_S_Se_ZVD < GP/C60_O_ZVD < GP/C60_Se_ZVD, except for GP/C60_S_ZVD, whose values increased when it was doped from 0.527 eV to 1.495 eV. Conversely, in terms of chemical hardness, an increase is noted after doping: the complex hardness increased except for that of GP/C60_S_ZVD, whose values decreased from 0.947 eV to 0.334 eV. This indicates an enhanced instability and reactivity of the complex after interaction with the ZVD, as evident in the values presented in the table. Furthermore, the ionization potential values determine the stability and reactivity, where the value of the electron affinity provides insight into the types of bonds and electron transfer, and the ionization potential signifies the interaction and stability [57]. In this study, the value of the IP during before interaction with the ZVD molecule shows an increasing trend of GP/C60_O < GP@C60_Se < GP@C60_S < GP/C60_O_S_Se. Notably, the value before interaction is less than the the complexes value, with the exception of GP/C60_S_ZVD, where the surface is higher than the complex by 5.278 eV, while GP/C60_O has the least IP before interaction by 5.174 eV. Upon complexing, the trend in its ascending order of GP/C60_S_ZVD < GP/C60_O_S_Se_ZVD < GP/C60_Se_ZVD < GP/C60_O_ZVD provides insight into the potency of the complex as a promising candidate for drug delivery, as GP/C60_S_ZVD possesses the least IP value, making it very reactive and easily releasing electrons for possible reactions in the environment. The electronegativity of GP/C60_S_ZVD, with a value of 20.209 eV, further explains why the complex possesses the properties needed for effective drug delivery [58]. It can be deduced from the results that there is a strong correlation between the calculated parameters and the delivery of zidovudine to the GP/C60_Se surface because of the greater energy gap upon interaction with the drug.

Fig. 3

Visual representation of the HOMO-LUMO isosurfaces of GP@C60_O_ZVD, GP@C60_O_S_Se_ZVD, GP@C60_S_ZVD, and GP@C60_Se_ZVD depicted in (a), (b), (c) and (d), respectively

Table 2 HOMO-LUMO results and quantum parameters such as ionization potential (IP), electron affinity (EA) chemical hardness (η), chemical softness (σ), chemical potential (µ), electrophilicity index (ω), and Fermi energy level (EFL)

Second-order perturbation energy analysis

Natural bond orbital (NBO) analysis provides a localized and chemically intuitive description of the molecular electronic structure, explaining intermolecular charge transfer [59]. Table S3 contains information on the transitions, donor and acceptor atoms, stabilization energy values E⁽²⁾, energy differences (E(j)-E(i)), and NBO overlap values (F(i, j)) for various transitions in different compounds (GP@C60, GP@C60_O, GP@C60_S, GP@C60_Se, GP@C60_O_ZVD, GP@C60_S_ZVD, GP@C60_Se_ZVD, GP@C60_O_S_Se_ZVD). This study explains the interaction between the compound and the doped compound as well as the interaction with ZVD. The transitions involving different orbital interactions, such as π→π*, π^*→π^*, LP→d*, describe the complicated nature of chemical bonding [60]. Table S3 shows the stabilization energy level obtained during the transition, where energy differences (E(j)-E(i)) explained energy changes associated with specific electronic rearrangements and NBO overlap values (F(i, j)) quantify the extent of orbital overlap, providing a crucial parameter for understanding the strength of chemical bonds in these compounds [61]. According to the study, the greater the value of the donor-acceptor interaction of the stabilization energy E⁽²⁾, the more stable the system. The GP@C60 compound demonstrates π→π* transitions involving various donor and acceptor atoms, revealing the involvement of specific carbon atoms in these bonding interactions with stabilization energies of 10.08 kcal/mol, 9.97 kcal/mol, and 15.02 kcal/mol. Additionally, the stabilization energies E⁽²⁾ of the interactions GP@C60_O, GP@C60_S, and GP@C60_Se were 98.21, 75.67, and 65.20, 80.64, 55.35, and 46.19, respectively, and the stabilization energies of the compound GP@C60_O were 98.21, 75.67, and 65.20, respectively. Moreover, the stabilization energies of the compound GP@C60_S were 80.64, 55.35, and 46.19, while those of GP@C60_Se was 84.50, 93.21, and 76.22, respectively. However, the complexes that interact with ZVD are GP@C60_O_ZVD, GP@C60_S_ZVD, GP@C60_Se_ZVD, and GP@C60_O_S_Se_ZVD, which are 43.74, 21.83, 60.01 and 96.44 kcal/mol; 75.87, 47.76 and 66.98, 96.17, 71.40 kcal/mol; and 80.38, 76.42, 65.42 kcal/mol, respectively. Upon interacting GP@C60 with oxygen, sulfur, and selenium, there was a notable increase in the stabilization energy of the system, indicating intense intermolecular charge transfer. Furthermore, it is worth noting that the system again interacted with ZDV, and there was both an increase and a decrease in the stabilization energy. Notably, the overall prominent transition that took place in the system was π$\:\to\:$π^* before the interaction with ZVD, and in the other way, the prominent transition changed from π$\:\to\:$π^* to π^*$\:\to\:$π^*. Hence, the major transitions occurred with π^*$\:\to\:$ π^*, contributing to the stability characteristic observed in the systems.

Density of states (DOS) analysis

The density of states (DOS) is employed to simulate and predict the electronic structure of materials accurately and determine the distribution of electrons among the individual contributing molecules [62]. Similarly, in quantum mechanics, the DOS represents the number of electronic states per unit energy interval at each energy level [63]. Different density of states parameters is employed in the determination of the electronic properties of the studied compound, which are the total density of states (TDOS), overlap partial density of states (OPDOS), and partial density of states (PDOS). The studied compound is composed of carbon (C), hydrogen (H), sulfur (S), oxygen (O), nitrogen (N), and selenium (Se) atoms, with each atom represented by a colour notation. From the results of this study, as shown in Figure S1, the carbon atom shows the highest orbital contribution in all the systems. Moreover, the Fermi energy levels of -6.78 eV, -6.78 eV, -6.89 eV, and − 6.89 eV for GP/C60_O_ZVD, GP/C60_S_ZVD, GP/C60_Se_ZVD, and GP/C60_O_Se_ZVD, respectively, are expressed with dotted lines. This shows non-dispersible changes in the energy fermi level of the compounds. The DOS is vital for designing new materials with specific electronic properties or understanding the behaviour of existing materials under different conditions because of the utilization of the studied nanocomposites in the delivery of drugs.

Visual studies

Analysis of the quantum theory of atoms-in-molecules (QTAIM)

QTAIM analysis was employed to thoroughly investigate the nature of the interactions between atoms within a molecule. This approach was applied to understand how interactive and stable our studied delivery material would be following interaction with the understudy ZVD [64, 65]. QTAIM studies have gained wide acceptance in drug delivery studies because they encompass important parameters that reveal different interactive behaviours via bond critical points [66]. The AIM molecular graphs of GP/C60_O_ZVD, GP/C60_S_ZVD, GP/C60_Se_ZVD, and GP/C60_O_Se_ZVD are shown in Figure S2, while the statuses of all the bond points and bond critical points related to these complexes are indicated in both the figure and Table 3. As presented in the table, all the complexes exhibited a large increase in the Laplacian charge density, demonstrating a decrease in the charge density between the adsorbent and the ZVD fragments. Thus, these trends present the typical characteristics of closed-shell interactions, indicating noncovalent interactions [67]. The value of ∇2 ρ(r) at the intermolecular hydrogen bond critical points of the GP@C60_S_ZVD complex was greater (0.0752 a.u. at H124-H131 bonding) than that of the other complexes, which confirmed that the greatest electronic change occurred in the nanoparticles despite intermolecular interactions with ZVD. This suggests enhanced stability and potentially stronger interactions, which are crucial for maintaining the integrity and efficacy of drug delivery systems. Conversely, the localization and delocalization of electrons within a molecule were investigated through the electronic localization function (ELF), where a value less than zero indicates high delocalization. In contrast, a value greater than zero indicates localization regions [68]. Based on our results, the overall complexes exhibit delocalization of electrons at most of their critical points. However, the lowest values were observed for S93-N107 and 217a.u. for GP@C60_O_S_Se_ZVD; for O96-N105 and 181a.u. for GP@C60_O_ZVD; for C18-N110 and 176a.u. for GP@C60_S_ZVD; and for C49-N120 and 131a.u., respectively. This implies that the electrons of atoms within these bonds are easily delocalized. Despite having ELF values < 0, the GP@C60_S_ZVD complex had higher values than the other complexes. Therefore, higher ELF values observed at most of the critical point implies strong interactions between the drug and carrier, promoting stability and controlled release dynamics. Additionally, the stability of the interactions can also be accounted for by the interaction between electrons and nuclei of the system, based on the values of potential energy density V(r), as higher negative values indicate significantly stronger interactions. Herein, the V(r) values were negative across all the critical points recorded for all the complexes. This suggests attractive interactions between the drug and the carrier, which can enhance stability by promoting complex formation.

Table 3 QTAIM results for the investigated complexes showing the critical point (CP), density of all electrons ρ(r), Laplacian of the charge density ∇2 ρ(r), kinetic energy density G(r), hamiltonian kinetic energy K(r), potential energy density V(r), energy density H(r), and electron localization function (ELF)

Noncovalent interaction studies

NCI serves as a critical tool for assessing the tenacity of intra- and intermolecular interactions, providing insight into the bonding between drug delivery systems and drugs, especially in the context of targeted therapies. However, in drug delivery studies, the primary goal is to transport therapeutic compounds to their intended site of action with precision and efficiency while minimizing off-target effects. Therefore, achieving this objective requires the evaluation of interactions between the drug and delivery material; hence, this study highlights the different natures of interactions/bonds that exist in the studied complexes, ranging from van der Waals interactions to hydrogen bond interactions and steric repulsion [69]. The study utilized a Multiwfn analyzer and visual molecular dynamic visualization tool to plot and extract the colours of the kinds of interactions present. The literature has reported that blue, green, and red spatial fragments represent different kinds of intermolecular interactions, where red signifies strong interactions solely due to steric repulsion, the green zone represents weak van der Waals forces within the studied compound, and the blue zone represents hydrogen bonding [70,71,72]. As noted from the results of the analysis for this investigation in Fig. 4, all four interactions exhibit large green isosurfaces within their intermolecular interactions. This reveals a weak force of interaction, highlighting the ability of the delivery material to easily release the ZVD compound. This effect was notable between the two combined nanoparticles, indicating that an electrostatic force existed between them, which could be the result of the properties of the atoms present and their conformational arrangements. Likewise, a large red band was observed for the intramolecular delivery system, while a few red signals were observed for the intramolecular delivery system of the ZXD compound, indicating the presence of strong and weak steric effects, respectively. Generally, this implies that the three-dimensional arrangement of atoms within a molecule significantly impacts its reactivity, stability, and behaviour due to the repulsive interactions between electrons and nearby atoms.

Fig. 4

Noncovalent interaction analysis identifying the nature of intermolecular interaction

Adsorption energy/BSSE analysis

Adsorption energy is an important aspect of the interaction between adsorbates and surfaces and plays an important role in understanding and optimizing processes in different fields, for instance, in the fields of catalysis gas storage and environmental remediation [73]. The adsorption energy is typically expressed in units of energy per absorbate molecule in kJ/mol or eV, as shown in Table 4, and the adsorption energy is also the energy change associated with the adsorption process. The adsorption energy was computed using the following equation:

$${{\rm{E}}_{{\rm{ads}}}}\,{\rm{ = }}\,{\rm{Ecomplex}}\,{\rm{ - }}\,\left( {{\rm{Eadsorbate}}\,{\rm{ + }}\,{\rm{Esurface}}} \right)$$

where E_ad is the adsorption energy, E_complex is the energy of the interaction (adsorbate on the surface), E _adsorbate is the energy of the isolated adsorbate and E_surface is the energy of the clean surface. Adsorption energy analysis can be said to be more favorable if the reaction process is exothermic (negative adsorption) < 0, while an endothermic process (positive adsorption) ≥ 0 implies less favourable adsorption [74]. Adsorption energy is a valuable tool that can be applied to predict the stability and efficiency of the adsorption process. A strong adsorption energy often corresponds to a stable configuration, resulting in a strong interaction between the adsorbate and the surface. However, weak adsorption, which in this instance is a positive adsorption energy value, may result in reversible adsorption, making surface properties important for a specific application. According to our results, all the complexes exhibit negative adsorption energies (GP/C60_S_ZVD − 1.599 eV, GP/C60_Se_ZVD − 0.400 eV, GP/C60_O_ZVD − 0.403 eV, and GP/C60_O_S_Se_ZVD − 0.390 eV), implying exothermic interactions. In other words, the complexity between the studied delivery materials and ZVD results in better adsorption, which is an important tool for drug delivery. Even at basis set superposition error (BSSE) evaluation, the systems remain consistent at the same adsorption phenomena, thereby demonstrating a tenacious adsorption capacity. Comparatively, the graphene-doped fullerene (C60) along with doped sulfur in the GP/C60_S (-1.599 eV) complex is considered suitable for drug delivery, especially for ZVD, and this more exothermic process is attributed to the conformation of the delivery system, where sulfur plays a significant role in surface engineering.

Table 4 The adsorption energy result, including E_complex (complex energy), E_adsorbate (adsorbate energy), E_adsorbent (adsorbent energy), given in Hartree units and Basis set superposition error (BSSE)

Drug release, pH, and dipole moment studies

The effectiveness and efficacy of how the nanomaterial will deliver the drug to the target site were determined by determining the parameters of the drug release mechanism, whereas the dipole moment was studied to determine the nature of the charge distribution within the system, and a pH evaluation showed how the environmental parameters influence the nanoparticle and how zidovudine (ZVD) can be released in the system [75,76,77]. This study investigated the effect of PH on zidovudine (ZVD) released from GP/C60_S, GP/C60_Se, GP/C60_O, and GP/C60_O_S_Se, which is a theoretical aspect of drug delivery in living systems. The impact of pH on drug release is a critical factor in drug development because it can significantly affect the solubility, stability, and release kinetics of pharmaceutical compounds. The study aimed to elucidate pH-dependent behaviours and their implications for drug release mechanisms from GP/C60-based systems. The pH influences drug release through its effect on the solubility and ionization state of both the drug and the carrier materials. At different pH levels, protonation or deprotonation of functional groups on the carrier surface can alter electrostatic interactions with the drug molecules, affecting their release kinetics [78]. Therefore, drug release values reported in atomic units (a.u.) provide insights into the release potential of zidovudine (ZVD) from different formulations. Interestingly, the drug release values varied across the various systems, as presented in Table 5, suggesting that the nature of the adsorption energy and the specific interactions between the drug and its carrier play a role in modulating drug release. Notably, GP/C60_S_ZVD exhibited the highest drug release value (-0.578 a.u.), indicating a tendency towards a less constrained release mechanism compared to GP/C60_Se_ZVD (-0.487 a.u.), GP/C60_O_ZVD (-0.515 a.u.), and GP/C60_O_S_Se_ZVD (-0.509 a.u.), suggesting potentially easier desorption under similar conditions. Furthermore, an examination of the dipole moment was conducted for each system to determine the nature of the charge distribution within the investigated systems and how drug molecules interact [79]. The dipole moment values were determined to be 6.480D, 4.203D, 3.804D, and 3.012D for GP/C60_S_ZVD, GP/C60_Se_ZVD, GP/C60_O_ZVD, and GP/C60_O_S_Se_ZVD, respectively. This suggests that the GP/C60_S_ZVD (6.480D) complex is more unlikely to separate due to the electrostatic forces between the adsorbate and adsorbent, as indicated by noncovalent interaction studies. While other systems with a dipole moment of less than 5D undergo quick desorption.

Solvation energy studies

In the development of a suitable drug delivery system, solvation energy explains the dynamics of solute-solvent interactions, which can be implemented in the field of drug delivery to assess the interaction between a drug and delivery material [80, 81]. The solvation energy denotes the energy variation that occurs when a drug molecule interacts with the molecule of delivery material. This energy change signifies the strength and nature of the bond formed between them. The solvation process typically involves the breaking of existing nanoparticle-nanoparticle and drug‒drug interactions, followed by the formation of new nanoparticle-drug bonds [82]. This rearrangement of molecules results in an energy change within the system. The solvation energy, therefore, quantifies this energy change, providing a measure of the affinity between the nanoparticle and drug [83, 84]. A solvation ≥ 0 energy indicates an endothermic process, meaning that energy is adsorbed as the nanoparticle interacts with the drug. This suggests a weaker interaction between the nanoparticles and drug molecules. Conversely, solvation < 0 energy denotes an exothermic process, where energy is released during solvation, implying a stronger attraction between the nanoparticles and drug molecules. The observed solvation energies in the systems are predominantly small and negative, indicating favourable interactions with the Zidovudine (ZVD) drug (Table 5). Notably, GP/C60_S_ZVD and GP/C60_O_ZVD exhibit the most negative solvation energies, measured at -0.825 eV and − 0.713 eV, respectively. This finding underscores the enhanced solvation interactions with ZVD brought about by incorporating sulfur (S) and oxygen (O) atoms in co-surface engineering, suggesting that these systems are highly useful for ZVD delivery. Conversely, GP/C60_Se_ZVD and GP/C60_O_S_Se_ZVD display marginally less negative solvation energies, recorded at -0.662 eV and − 0.681 eV, respectively. Although still indicative of favourable interactions, these values suggest that the initial GP/C60_Se_ZVD system and its derivative with added S and O atoms may exhibit slightly weaker solvation interactions with ZVD than systems lacking selenium (Se) atoms. This phenomenon might arise from the broader electronegativity of selenium, contributing to the relatively less attractive characteristics of the investigated delivery systems. In general, the incorporation of specific atoms, such as S, O, and Se, in the co-surface engineering of graphene-based systems has visible effects on solvation energies, influencing the attractiveness and effectiveness of these systems for ZVD delivery. Understanding these dynamics aids in the rational design and optimization of custom-made drug delivery platforms for specific therapeutic agents.

Table 5 Results from the analysis of drug release (a.u.), solvation energy (eV), and dipole moment (D)

Conclusions

Encapsulating antiretroviral drugs within nanoparticles is an excellent innovation for precise drug delivery. Therefore, in this study, we re-engineered nanoparticles by combining two surfaces, namely, graphene and fullerene C60, and modified with sulfur (S), selenium (Se), and oxygen (O) for the delivery of zidovudine (ZVD) using density functional theory (DFT) simulations, which were calculated with the B3LYP functional and D3(BJ)/Def2svp basis sets. Moreover, the modification of the aforementioned elements could enhance the attractiveness of the studied delivery material, thereby facilitating the delivery of ZVD. In our research, a structural study was performed to determine the structural changes that occurs on customized delivery materials and their interactions with ZVD. Most of the selected bond labels displayed decreased bond lengths after interaction with ZVD, indicating stable interactions. However, the after-bond length of GP@C60_S_ZVD was shorter than that of the other materials, implying better stability. Further evidence of a stable nature was obtained in second-order perturbation energy analysis where the π^*$\:\to\:$π^* transition is predominant among the systems. In addition, the GP/C60_S_ZVD (0.668 eV) and GP/C₆₀_O_S_Se_ZVD systems exhibited favourable interactions following the HOMO-LUMO energy gap. The energy gap showed decrease in the order GP/C60_Se_ZVD (2.083 eV) > GP/C60_O_ZVD (1.778 eV) > GP/C60_O_S_Se_ZVD (1.303 eV) > GP/C60_S_ZVD (0.668 eV). Visual studies, especially quantum theory of atoms-in-molecules (QTAIM) and noncovalent interactions (NCIs), have revealed the intramolecular and intermolecular interactions of the complexes. The findings show that all the studied systems involve noncovalent interactions rather than covalent interactions, providing good characteristics for ZVD delivery. Further investigations in NCI studies revealed intermolecular steric repulsion via intramolecular and van der Waals interactions across all the tested systems. The sensitivity of the delivery materials to ZVD was revealed through the Fermi energy density of state results. The Fermi energy levels of GP/C60_O_ZVD, GP/C60_S_ZVD, GP/C60_Se_ZVD, and GP/C60_O_Se_ZVD are − 6.78 eV, -6.78 eV, -6.89 eV, and − 6.89 eV, respectively. This shows a similarity in sensitivity capacity; hence, despite the variation in doped atoms, they possessed a similar nature in terms of sensitivity. Their level of interaction was further studied via adsorption energy, whereby GP/C60_S_ZVD, GP/C60_O_ZVD, GP/C60_Se_ZVD, and GP/C60_O_S_Se_ZVD had energies of -1.59949 eV, -0.40265 eV, -0.40007 eV, and − 0.39032 eV, respectively. The mechanism of drug release was investigated through pH and dipole moment, where GP/C60_S_ZVD was indicated to have the highest dipole moment (6.480074 D), suggesting quick release. However, the dipole moments of GP/C60_Se_ZVD, GP/C60_O_ZVD, and GP/C60_O_S_Se_ZVD are 4.203098 D, 3.804473 D, and 3.019110 D, respectively, which provide a significant model for the subsequent delivery of the zidovudine drug in cases of HIV/AIDs infection.