Introduction

The COVID-19 pandemic brought the need for advances in science, and manipulation of materials at the nanometer scale is highlighted, as it reveals extraordinary physiological and chemical properties [1]. Vaccines have been developed to prevent severe cases of the disease; however, non-pharmacological measures are still needed to slow down the spreadness and reduce morbidity and mortality from this disease [2]. In this scenario, drug nanocarrier systems have gained prominence in research, as the search for greater efficiency in the treatment of the virus has become essential. Thus, it seeks to avoid symptoms, serious consequences, and even death [1]. It allows effective and specific drug delivery to the target cells, which reduces the necessary dose and the possible side effects for the individual [3].

The SARS-CoV-2 virus causes this disease, which has in its composition a protein named spike, which gives it the classification of coronavirus. This protein is responsible for the cell internalization of the virus, as it is recognized by the human angiotensin-converting enzyme-2 (ACE-2) and within the cell undergoes replication controlled by the main protease 3-chymotrypsin-like cysteine (Mpro or 3CLpro) [1, 3]. This virus undergoes mutations that classify it into different variants that can cause different symptoms. Among the variants, the Delta variant has the highest affinity for interaction with ACE-2, and the Omicron variant has mutations that make it more likely not to be tackled by vaccines already in development [4]. Because of this, several studies are being carried out via in silico methodologies for screening compounds for interacting with the spike protein, ACE-2, and Mpro to develop alternative methods of diagnosis and treatment. In this way, Eskandari [5] carried out molecular docking studies to identify natural compounds as potential spike protein receptor-binding domain (RBD) and with Mpro inhibitors. The author suggested that vitamins such as riboflavin are potential inhibitors of Mpro, vitamin B12, and benfotiamine are inhibitors of RBD, and bentiamine and folic acid are potential inhibitors for both. In this line of research, the natural compounds such as kadsurenin L and methysticin [6] and p-coumaric acid, allagic acid, kaempferol and quercetin [7] are identified as Mpro inhibitors.

Additionally, carbon nanomaterials such as graphene, and their derivatives like graphene oxide (GO) and carbon nanotubes, which are composed of a honeycomb packed on sp2 hybridization, are widely used for nanomedicine applications [1, 8, 9]. These nanomaterials are already studied as nanocarriers for antiviral drugs against HIV and reovirus [10]. GO is derived from graphene, with the presence of functional groups that facilitate its functionalization, that is, the interaction with other molecules of interest. It has an atomic thickness, which gives it the 2D structure classification and features such as biocompatibility, energy transfer efficiency, and amphiphilicity [11]. These properties are important for biomedical applications of imaging, biosensors, and drug nanocarriers [12]. Also, GO has antiviral activity by itself and demonstrably inhibits viral infection of a type of coronavirus that infects pigs through structural destruction [13]. Several studies also show that graphene and graphene oxide are more efficient than other nanomaterials [14], with regard to reducing the transmissibility and infectivity of SARS-CoV-2, being promising for use in prophylactic approaches, such as masks or surface coatings [14, 15].

Other nanomaterials such as molybdenum disulfide (MoS2) are promising alternative therapeutics to inactivate SARS-CoV-2. Bisht and collaborators [16] evaluated the interaction of 2D MoS2 nanosheets with the SARS-CoV-2 spike protein, the human ACE-2 receptor, and the complex formed between them through molecular docking and atomistic simulations. The results show that MoS2 nanosheets can effectively bind to the receptor-binding domain (RBD) of the spike protein with good docking energies, and the interactions occur through various hydrogen bonds and van der Waals interactions. Thus, the authors clearly demonstrate the antiviral potential of MoS2 2D nanosheets [16].

In parallel, flavonoids are natural compounds with excellent pharmacokinetic activity, which have anti-inflammatory, antioxidant, bactericidal, and antiviral activities, and the ability to boost the immune system. Because of this, they are presented as proposals in the fight against COVID-19, given the urgent need to remodel drugs [17, 18]. Lakshmi and collaborators [19] point to orientin as one of the natural compounds that are presented as a proposal to alleviate the viral infection, without side effects, by promoting the immune system. His research involved in silico investigation via molecular docking with the spike protein and the main protease of the virus, as well as with the human ACE-2 receptor. Additionally, Gentile and collaborators [20] investigated molecular dynamics, where they point to apigenin as one of the compounds capable of inhibiting the main protease of SARS-CoV-2. About nanocarriers, Rahmanian and collaborators [21] synthesized nano GO as a carrier of quercetin, an antioxidant flavonoid capable of promoting the immune system, and point out that the nanomaterial is an effective, low-cost option with high production capacity for interaction with this type of biomolecule, especially for oral internalization.

In this paper, first-principles computer simulations were performed to study the interaction between graphene oxide and the biomolecules orientin and apigenin. It is intended to analyze the stability of the proposed system to link the therapeutic activities of these flavonoids to a nanocarrier system as a treatment method against COVID-19. After that, molecular docking studies were carried out to identify the interaction of the graphene oxide, apigenin, and orientin structures with the RBD of the spike protein of two coronavirus variants, Delta, and Omicron.

Methods

Firstly, the interaction between GO and the biomolecules apigenin and orientin was studied through first-principles simulations based on Density Functional Theory (DFT) [22] to evaluate the proposed structures for a nanocarrier system against COVID-19 by their electronic properties. Then, a molecular docking study was performed with a target in the RBD of the spike protein of the Delta B.1.617.2 and Omicron B.1.1.529 variant to identify the conformation and the protein–ligand affinity [23]. The RBD of the spike protein was specifically chosen for its importance in the recognition of the protein by the cell and consequently in the viral infection.

DFT simulation

DFT simulation and calculation parameters

The ab initio methodology allows the understanding of chemical, physical, and biological phenomena by solving Schrödinger’s equation, which describes the electron movement around the atom nucleus by their wave function. Meanwhile, the DFT is a theory formulated by Hohenberg and Kohn, which reduces the computational effort by indicating theorems that show that the electron density can take the focus instead of each electron of the wave function.

The computational simulation was performed using the SIESTA software (Spanish Initiative for Electronic Simulations with Thousands of Atoms) [24], which uses the DFT [25] and some approximations to solve the electronic problem. The local density approximation [26] was used as an auxiliary tool for computational calculations, as it provides us with theoretical results equivalent to the experimental ones in weakly interacting systems, such as simulations involving GO and organic molecules [27]. Pseudopotentials were used to simplify the description of electrons. The double zeta polarization (DZP) basis set was chosen, with an energy shift of 0.05 eV and energy cutoff of 200Ry, and all the calculations were performed to converge when the residual forces on each atom were smaller than 0.05 eV/Å, for all the simulations. The super-cell size was 30 Å in x, y, and z directions, which is enough size so that the replicated system does not interfere with the simulation; that is, it does not become periodic. The methodology is similarly approached on other studies, which verified the adsorption of small aromatic molecules in pristine and functionalized graphene [28], of 17 β-estradiol in graphene oxide [27], and of pentachlorophenol on carbon nanostructures [29], and also, to identify a weak interaction force between graphene and doxorubicin [30], and the stability of π-aromatic interaction of different groups adsorbed in graphene surface [31].

The isolated biomolecules apigenin (PDB: 5280443) and orientin (PDB: 5281675), and GO interacting with each of the biomolecules, were analyzed separately. We have analyzed the gap between HOMO (highest occupied molecular orbital) and LUMO (lowest unoccupied molecular orbital) of the system. The charge plots of them were visualized by choosing 0.001 eV/Å3 as the contour value. The binding energy was calculated according to Eq. 1 which calculates a basis set superposition energy (BSSE) [32]. By BSSE, the interaction energy is calculated considering the initial configuration of the formed interacting system (AB). For this, it is subtracted from the energy by A without explicitly considering B, which is treated as a ghost (A + B(ghost)), and vice versa (A(ghost) + B) [31].

$$E_{\mathrm{binding}}\;=\;E_{AB}\;-\;E_{A(\mathrm{ghost})+B}\;-\;E_{A+B(\mathrm{ghost})}$$
(1)

The charge transferences between all systems are calculated as shown in Eq. 2 for the biomolecule, so the negative/positive transferences indicate that the biomolecule donates/accepts electronic charge.

$$\mathrm{Charge}\;\mathrm{transference}\;=\;E_{\mathrm{theoretical}}\;-\;E_{\mathrm{resultant}}$$
(2)

These analyses will show us through the energy and energy difference in eV, and the characteristics of interaction between GO and biomolecules. Thus, it is possible to assume that the proposed system fulfills its characteristics as a potential nanocarrier.

Molecular docking

Structures preparation, molecular docking, and results analysis

The SARS-CoV-2 Delta variant Crystal structure file (PDB:7V7V) and the SARS-CoV-2 S Omicron Spike (PDB:7QO7) were treated using the AutoDockTools software, which required the addition of charges and hydrogen atoms, where the conversion to pdbqt format was also performed, while the three-dimensional structures of the ligands apigenin and orientin were acquired from PubChem, treated by adding a torsion point in the structure, and converted to pdbqt format also by AutoDockTools [33].

The molecular docking was performed using AMDock software, which is a graphical tool that helps docking between protein and ligand using AutoDock4 [34]. Autodock proved to be an effective tool in predicting binding conformation and binding energy of flexible ligands to flexible macromolecules through a semi-empirical binding affinity force field. Autodock 4 software uses a scoring function to evaluate interactions, returning us a binding affinity value in kcal/mol, by which it is possible to evaluate potential inhibitors and binding states [35]. The ligand-docking calculation follows Eq. 3.

$$\mathit\triangle G\mathit=\mathit{\left({V_{bound}^{L-L}-V_{unbound}^{L-L}}\right)}\mathit+\mathit{\left({V_{bound}^{SP-SP}-V_{unbound}^{SP-SP}}\right)}\mathit+\mathit{\left({V_{bound}^{L-SP}-V_{unbound}^{L-SP}}\right)}\mathit+\mathit\triangle S_{\mathit c\mathit o\mathit n\mathit f}$$
(3)

where L refers to the “ligand,” SP refers to the “spike protein,” ∆Sconf is an estimate of the conformational entropy lost upon binding, and V is a pair-wise evaluation for dispersion/repulsion energies, hydrogen bonds, electrostatics, and desolvation [36]. These binding affinity analyses are expected to be useful in rational drug design by virtual screening and other applicability. The default parameters for the grid box were centered at 221.3 Å × 175.5 Å × 276.3 Å (x, y, and z) to cover the region of RBD of the Delta spike protein [37]. For the Omicron spike protein, the default parameters for the grid box were centered at 213.6 Å × 181.5 Å × 134.3 Å (x, y, and z) to cover the region of RBD (residues 330–520) [38].

The results were analyzed regarding the conformation score, the interaction affinity resulting from the best conformation, that is, the one with the most negative binding energy, and the types of protein–ligand interactions formed. The stability of the docked structures was examined using the RMSD (root-mean-square deviation) parameter, which is the distance between the ligand and receptor and must be less than 2 angstroms [33]. Discovery Studio 2021 software was used to visualize the interaction results and the 2D map of interaction was performed through LigPlot software. The proposed methodology was similarly approached in other studies, such as the interaction in the complex formed by graphene and glutamate [39], and for single-walled carbon nanotubes and human mitochondrial voltage-dependent anion channel [40].

Results and discussion

This study aimed to analyze, by means of first-principles calculations based on DFT, the feasibility of interaction between the GO, and the flavonoids apigenin and orientin, for the development of a nanocarrier; and, by means of molecular docking, the interaction of these structures to the spike proteins of the Delta and Omicron variants separately, to identify the ligands as inhibitors of the target protein, therefore, acting against the SARS-CoV-2 viral infection. All results showed RMSD lower than 2 angstroms; therefore, they are considered reliable docking results. In summary, the stability of the system formed with GO as a nanocarrier of the biomolecules apigenin and orientin was guaranteed by ab initio simulation, while the interaction of ligands with target proteins was verified by molecular docking.

DFT simulations and analysis

Isolated structures

Initially, we evaluated the isolated structures of GO, apigenin, and orientin (Fig. 1). The energy levels/bands and the plot of the charge in the region of HOMO (highest occupied molecular orbital) and LUMO (lowest unoccupied molecular orbital) were analyzed. The studied GO has the molecular formula C55H21O6; it is a planar structure of graphene with chemical groups functionalizing the structure, such as carboxyl, hydroxyl, and epoxy, which allows the structure to be more interactive and provides for biological applications [11]. The results show that the GO has a HOMO/LUMO difference of 0.60 eV, similar to the results of De Oliveira and collaborators [27]. In this case, the charge plot on HOMO is arranged at the ends of the sheet, with emphasis on the carboxylic group, and on LUMO in the central region, especially on the hydroxyl groups, which shows us the greater interaction propensity of oxygen atoms. Furthermore, studies indicate that GO can be doped with aromatic molecules by interacting in regions containing COOH groups and via π-π stacking, taking advantage of their sp2 hybridization sites [41].

Fig. 1
figure 1

Structures studied (a). Energy bands and charge plot (b) in the HOMO and LUMO of the studied structures: GO, apigenin, and orientin. Value used for the isosurface is 0.001 e−/A3

Full size image

In the case of apigenin, it has a HOMO/LUMO difference of 2.50 eV, a value that agrees with the flavonoid gap variation result from Bitew and collaborators [42] who studied it by DFT/B3LYP hybrid functional. Additionally, the charge plot in HOMO is arranged over the aromatic part, while in LUMO it is arranged over most atoms, as is also found in the work of Zheng and collaborators [43]. On the other hand, orientin has a HOMO/LUMO difference of 2.50 eV, like the results of Li and collaborators [44], who used the Becke exchange methodology plus Lee-Yan-Parr correlation (BLYP). The charge plot of orientin is found in the literature, distributed over the aromatic part of the molecule in both HOMO and LUMO.

After evaluating the structural and electronic properties of the isolated structures, four interaction configurations were selected. Initially, all configurations were performed to evaluate the most favorable distance for the interaction to occur, with the distance up to 3 Å (ranging from 0.25 to 0.25 Å). When the best distance was found, the system was optimized for further analysis of the results. In the next topics, we will present the results of the interaction of GO with apigenin and with orientin.

GO and biomolecule interaction—DFT simulations and analysis

For the study of GO interacting with apigenin, we evaluated four distinct configurations (Fig. 2), where four configurations were assembled with the biomolecule perpendicular to the GO on the edge of the structure (GO + apigenin_1) and in the middle of the GO (GO + apigenin_2) and with the two structures parallel to each other, on the edges of the GO (GO + apigenin_3 and GO + apigenin_4). For the study of GO interacting with orientin, the configurations studied are shown in Fig. 3; the configurations were modeled so that the charge density region of the biomolecule remained parallel to the GO, approaching a region of few functional groups (GO + orientin_1), approaching the oxygen-containing groups of both structures (GO + orientin_2 and GO + orientin_3), and approaching the carbons of the aromatic rings of both structures (GO + orientin_4). Figure 3 shows the configurations studied for the interaction of GO with orientin. Table 1 shows the values of connection energy, distance, HOMO/LUMO difference, and electronic charge transfer obtained in the different configurations. The most stable results are marked, that is, the ones with the lowest energy. It is noteworthy that all the resulting energies are negative and, therefore, indicate an interaction between the systems.

Fig. 2
figure 2

Studied configurations for GO and apigenin interaction (a). Energy bands of GO, apigenin, and GO + apigenin_4 (b). Charge plot in HOMO and LUMO for the most stable configuration (value used for the isosurface is 0.001 e−/Å3) (c)

Full size image
Fig. 3
figure 3

Studied configurations for GO and orientin interaction (a). Energy bands of GO, orientin, and GO + orientin_3 (b). Charge plot in HOMO and LUMO for the most stable configuration (value used for the isosurface is 0.001 e−/A3) (c)

Full size image
Table 1 Results of the interaction in different configurations of the system containing GO and apigenin
Full size table

The simulations demonstrate that the perpendicular configurations resulted in the least attractive energies, while the system GO + apigenin_4 is the configuration of greater stability because it comparatively presents a lower value of total energy. In this configuration, apigenin is arranged in a planar way over the GO, in which we obtained binding energy of −0.96 eV. In this system, the HOMO/LUMO difference was 0.53 eV, like GO alone. Thus, it is also verified that in the energy bands of GO + apigenin_4, which are shown in Fig. 2b, the association of apigenin does not cause considerable modification in the bands of isolated GO.

The calculations of the most stable distances identified that the most stable system was the GO + orientin_3, as it comparatively presents a lower total energy value. This configuration resulted in the binding energy of −1.86 eV. The energy bands are shown in Fig. 3b, where we can see that the GO + orientin_3 system presents a HOMO/LUMO difference of 0.36 eV, a value that shows some modification when compared to the isolated GO (0.61 eV). Additionally, coupled with the low value of total energy, it indicates a physical adsorption system, where the association of compounds does not cause chemical modification of them; however, interaction still occurs due to negative energy values.

The least attractive configurations which occur in configurations 1 and 2 can be explained by the repulsions between the hydrogen atoms [31]. A theoretical study [45] of the interaction of apigenin with hexagonal boron nitride (h-BN), a structure similar to pure graphene, through DFT with the Lee–Yang–Parr [46] functional through the Dmol 3 [47] (B3LYP) package shows that the interaction is favored by “stacking-AA.” In the cited work, the authors obtained binding energy in the order of 1.61 eV, like that found in this work, of 1.24 eV, emphasizing that they are different structures, only similar. The author suggests that the interaction is favored by the interactions of the aromatic rings of apigenin through hydrophobic interactions and π-π bonding interactions.

The calculated charge transfer was 0.21 e− for apigenin, which acts as a charge donor. The HOMO and LUMO charge density plot is shown in Fig. 3c, where we observe a concentration of electronic charges for both HOMO and LUMO only on GO and the absence of charges on apigenin. The interaction of GO occurs by physical adsorption with apigenin, and this is a weak interaction. The calculated charge transfer was 0.03 eV; that is, orientin receives charge from the GO. The loading plot for the GO + orientin_3 configuration shows that both in HOMO and LUMO the electronic charge is localized on the GO as it is shown in Fig. 3c, indicating that a weak interaction occurs, that is, physical adsorption. We emphasize that the simulations studied in both systems showed better results when configured by keeping the GO and the biomolecule in a parallel position with each other, which is due to attractive π-π staking [31, 41].

Molecular docking simulations and interaction analysis

The molecular docking results showed different conformations in relation to the target protein; however, the dockings with the best score were chosen for comparative presentation in Table 2. The 2D maps of interactions made by LigPlot are shown in Fig. 4.

Table 2 Molecular docking results of ligand interactions against RBD Spike protein of Delta and Omicron variants of coronavirus
Full size table
Fig. 4
figure 4

Binding sites of Delta (a) and Omicron (b) spike proteins, and two-dimensional map of interactions between the protein target and the proposed ligands (c)

Full size image

The analysis of GO exhibits four hydrogen bonds (H bond) in chain A (Tyr394, Arg464, Arg353 and Arg355) for Delta, and just one H bond (Ala519) in chain A for Omicron. Also, for the Omicron variant these interactions are made all in chain A, while GO interacts with Delta in chains A and B. Apigenin forms just one H bond with Delta chain A (Gln472), and four H bonds with Omicron in chains A (Tyr418 and Arg454) and C (Phe372 and Phe374). And orientin forms four H bonds, two of them on amino acid residue Leu333 chain A, and the other two on amino acid residue Glu167 chain B for the Delta variant, and for the Omicron variant, three H bonds on chain A (Thr330, Asn357, and Pro327). The non-ligand hydrophobic bonds (HP-bond) are all different when comparing the results between variants, so we can assume that they are interacting with different areas.

For both the Delta variant and the Omicron variant, GO was presented with the lowest affinity (−11.88 and −13.20 kcal/mol, respectively). In the literature, GO is identified as a structure with potential for interacting with the spike protein, with values from −10.5 to −8.3 kcal/mol, with docking related to another variation of the spike protein [48].

Regarding the biomolecules, both apigenin and orientin show themselves as better inhibitors against the Omicron variant of SARS-CoV-2 spike protein, presenting a significant number of hydrogen bonds and affinity binding of −9.6 kcal/mol and −9.58 kcal/mol, respectively. This must occur because the Omicron variant presents several modifications on its structure including on RBD which increases its transmission capacity [49]. The binding affinity reached for apigenin is in accordance with the studies of Subbaiyan (2020) [50] and also with the values reached for polyphenols catechin and curcumin [51] and for other flavonoids, which presented the affinity binding values from −11.2 to −8.0 kcal/mol against spike protein [52].

The binding affinity values are in agreement with results that indicate inhibitory activity of the spike protein with other binding structures [5, 48, 52,53,54 ]. From the point of view of binding energy, GO shows strong interactions with the targets of the variant Delta and Omicron of SARS-CoV-2.

The molecular docking points out to the structures GO, apigenin, and orientin as potential inhibitors of the RBD of spike protein of the Delta and Omicron variants presenting results similar to those pointed out by other authors [5, 48, 52,53,54]. These results coupled with other characteristics such as the large surface area of the GO allow it to act as a nanocarrier [9]. Studies show that the proposal of a GO-based nanocarrier doped with molecules via π-π bonding, as pointed out in this work, is successfully carried out [41]. The therapeutic activities of the chosen biomolecules are pointed out in the literature; however, it is noteworthy that new molecular docking studies [55,56,57] and subsequent in vitro tests should be performed to analyze the interaction of the system with the target proteins of SARS-CoV-2, as a proposed treatment against COVID-19.

Conclusion

In this paper, we evaluated the interaction of GO with the biomolecules of apigenin and orientin through first-principles calculations based on Density Functional Theory. The results show that the most stable systems of the interaction of GO with apigenin and orientin present binding energy of −0.96 and −1.86 eV, respectively. These more stable configurations show that the interaction between GO and the biomolecules occurs when the structures are arranged parallel to each other; that is, it occurs via π-π stacking found in different studies. Together, the results of electronic charge transfer and 3D plot in HOMO and LUMO corroborate that the weak interaction occurs through physical adsorption. The resulting binding energies indicate physical adsorption between the systems, which is desired in a nanocarrier complex in order to remain the structure of the adsorbed biomolecule. At the same time, we studied the interaction of the ligands GO, apigenin, and orientin in relation to the RBD of two variants (Delta and Omicron) of the spike protein of SARS-CoV-2, the virus that causes COVID-19. It was verified that the protein–ligand interaction suggested the structures as potential inhibitors by forming bonds with amino acid residues on the RBD. Among the studied structures, the GO presented the lowest binding affinity of −13.20 and −11.88 kcal/mol. Regarding the biomolecules, they presented values lower than −6.0 reaching to −9.6 kcal/mol and are also suggestive to the spike protein inhibition because these values are like the ones reached in the cited studies. It is noteworthy that the best results in binding affinity of the structures were in relation to the Omicron variant of the coronavirus, which has been identified as a variant in 2021. The molecular docking results show that all the studied structures are presented as potential inhibitors of SARS-CoV-2 Delta and Omicron infection, while the ab initio results show that the GO-based system is promising as a nanocarrier of the biomolecules apigenin and orientin.