Introduction

In recent decades, viral diseases including H1N1 influenza, SARS-COV-1, COVID-19, and Ebola have threatened the health of human society. In 2019, a new coronavirus disease called COVID-19 emerged in Wuhan, China and quickly spread worldwide, resulting in millions of casualties [1]. The disease is caused by SARS-CoV-2, which shares similarities with the SARS-CoV-1 virus and has a vital protease, 3CL^pro or M^pro, that plays a significant role in virus replication [1-4]. Therefore, this protease is a potential target for developing antiviral drugs against SARS-CoV-2 [4-5,6-9].

Despite efforts by pharmaceutical companies and research groups, only a few monoclonal antibodies, such as bamlanivimab–etesevimab, casirivimab–imdevimab, and sotrovimab, have been authorized by the US FDA under emergency use authorization for treating mild to moderate COVID-19 in high-risk patients [10,11].

Remdesivir is another drug authorized for hospitalized patients, but it requires administration by injection or infusion in a therapeutic setting [12]. In December 2021, two new oral antiviral agents, molnupiravir and nirmatrelvir/ritonavir (Paxlovid), gained emergency use authorization from the FDA for adult patients with mild to moderate COVID-19 who were at high risk of progression to severe COVID-19 [12-14]. Tocilizumab and baricitinib have also gained EAU from the FDA for treating hospitalized COVID-19 patients aged 18 or above [5,14-15]. Moreover, computer-based drug design is a promising approach to enhance new antiviral components using repurposed FDA-approved drugs or natural products. The SARS-CoV-2 main protease (M^pro) plays a vital role in viral replication and is a good target for therapeutic repositioning [6].

Honeybee natural products are rich in flavonoids/non-flavonoids and have been used to treat various diseases due to their low toxicity and high biological potency. Some of these natural products, such as 3,4,5-Tricaffeoylquinic acid, Kaempferol-3-O-glycoside, and (+)-Pinoresinol, have shown antiviral activity against various viruses [4]. In this study, the antiviral effect of 102 flavonoids/non-flavonoids of honeybee propolis natural products was investigated against the active site of SARS-CoV-2 M^pro using molecular docking. The effects of these compounds on the structure of M^pro were also investigated by 150 ns molecular dynamics (MD) simulations, and pharmacokinetic properties, toxicity, and drug likeliness of the selected compounds with the lowest binding energies were considered.

Experimental

Preparation of ligands and receptor

The period of the present research has taken approximately eight months. To prepare for the study, the 3D structures of the compounds were gained from databases such as Zinc 15, PubChem, and ChemSpider in SDF format, and then they converted these structures into PDB format using open babel [16]. In addition, the 3D crystal structure of the SARS-CoV-2 main protease (PDB code: 6LU7) was obtained from the Protein Data Bank [17]. The Autodock Tools (ADT) software (version 1.5.6) was used for pre-processing steps [18] including adding the polar hydrogens, atom types of all protein and ligands, and gasteiger charges. For the ligands preparation, the researchers set the torsions to be active to develop searching phase space. Finally, the ligands were saved in pdb format for AutoDock Vina.

Molecular docking via Autodock vina

The process of molecular docking was executed through the utilization of Vina 1.1.2 Windows 10 platform (64-bit) with Lenovo IdeaPad L340-15IWL (CORE i7, system memory: 256 GB SSD, 12 GB RAM). The active site of the SARS-CoV-2 Main Protease was subjected to a grid with a dimension setting, wherein a grid spacing of 1. 0 angstrom was implemented [19]. The remaining parameters were maintained in their default setting. Based on the outcomes obtained from Autodock Vina, it can be surmised that the compounds exhibiting the lowest binding energies were distinguished as the foremost potential inhibitors of the main protease of COVID-19. The visualization of ligands’ interactions with the target was successfully obtained using Discovery Studio v 4. 5.

Molecular dynamic simulation

The molecular dynamics (MD) simulation was applied to study the dynamic behavior of bonded ligands to receptors at atomic scales. Three chosen protein-ligand systems (from the second stage of docking study) were submitted to MD simulations through Gromacs package (version 2021) to investigate further interactions between ligands and the targets. In terms of the topology data, these information were gained using gromos53a6 force field for receptors and PRODRG server for ligands, respectively [20,21].

Upon solvation of the simulation boxes utilizing SPC/E water and subsequent addition of adequate quantity of counter ions, neutralization of all protein-ligand systems was achieved [21]. Subsequently, the energy optimization procedure for the boxes was executed employing the steepest descent algorithm until obtaining an energy of less than 10 kj/mol [22].

Furthermore, the application of periodic boundary conditions was adopted in all boxes and along all directions (X, Y, and Z). The NVT and NPT protocols were employed with temperature and pressure coupling conditions established at 310 K and 1 bar, respectively. The equilibration stage involved the utilization of the V-rescale thermostat and a parinello-rahman barostat for the purpose of temperature and pressure equilibrations, respectively. Van der Waals and electrostatic interactions were estimated within a cut-off distance of 1 nm. The LINCS algorithm was implemented to enforce the bond constraint among heavy atoms [23]. The leap frog algorithm [24] was employed to conduct 150 ns molecular dynamics simulations. The investigation of interaction energies between ligands and targets was carried out utilizing the Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) method [25].

The 2D molecular interactions, coupled with graphical depictions, were derived from the software tools Ligplot [26] and VMD (https://www. ksuiucedu/Research/vmd). Subsequently, the potential toxicity of the fungal metabolites was evaluated using the ProTox-II webserver [26]. Pharmacokinetics properties were also achieved by SwissADME server (http://www.swissadme.ch).

Results and discussion

Docking validation

To validate the docking process, the re-docking experiment of the N3 inhibitor onto the binding site of the COVID-19 M^pro (PDB code: 6LU7) was conducted. The results showed a favorable binding energy of -13.4 kcal/mol and the inhibitor was found to make 13 hydrogen bonds with the receptor (Figure 2). This suggests that the docking process was accurate and reliable, and provides confidence in the subsequent analysis of the natural compounds.

Molecular docking results

In this study, natural compounds derived from propolis as inhibitors have been chosen to test their antiviral potency against the main protease of SARS-CoV-2 via AutoDock Vina. Armed with the Docking results, it can be mentioned that six out of all natural compounds delineated the lowest binding energies to the viral receptor, including 6-Cinnamylchrysin, (+)-Pinoresinol, Culifolin, 5,4'-Dihydroxy-3,3'-dimethoxy-2-prenyl-E-stilbene, Teferin, and 2-Acetoxy-6-p-hydroxybenzoyl jaeschkeanadiol. The binding energies along with the 2D and 3D visualizations of these natural inhibitors’ interactions with the viral protease have been illustrated in Table 1 and Figures 2 and 3.

As presented in Table 1, (+)-Pinoresinol and 6-Cinnamylchrysin have bound to the active site of the viral receptor with the lowest binding energies of -10.1 and -9.8 kcal/mol, respectively, which confirm their high antiviral potency. (+)-Pinoresinol has formed various interactions with the residues located in the active site, for example, the aromatic rings of this compound have formed Amide-π stacked and π-donor hydrogen bond interactions with the amino acids Gln189 and Cys145, respectively.

In addition, the compound has formed three hydrogen bonds with amino acids Leu141, Ser144, and Gly143, along with Van der Waals interactions, have stabilized the ligand conformation at the enzyme's active site.

Furthermore, these interactions can block Cys145, which is a nucleophile in the polyprotein cleavage reaction [27], and can also hinder the enzyme's binding to the substrate and prevent its proteolytic activity in cleaving polyproteins [27,7]. Based on docking data, it can be inferred that (+)-Pinoresinol can be a good inhibitor for blocking the viral enzyme (Figure 3A). Moreover, the docking results of 6-Cinnamylchrysin at the active site of the M^pro enzyme have been analyzed and shown in Figure 1B.

6-Cinnamylchrysin has been positioned at the active site of the protein by a H-bond with the His164 and a π-sulfur bond with His165, this compound may disrupt the formation of the enzyme's active site, as these two amino acids (in domain II) along with the residues Leu24, Pro39, and Met49 (in domain I) have vital responsibilities in generating the protease enzyme's active site [7]. In addition, the creation of π-alkyl and hydrogen interactions with the amino acid Cys145 may enhance the inhibitory potency of this compound against the protease enzyme.

   According to Figure 3C, D, it is worth mentioning that the other two natural compounds 2-Acetoxy-6-p-hydroxybenzoyl jaeschkeanadiol and 5,4'-Dihydroxy-3,3'-dimethoxy-2-prenyl-E-stilbene have shown the binding energies of -8.9 and -8.6 kcal/mol, respectively, showing their convincing antiviral effects against the viral enzyme. It can be mentioned that these two natural ligands have been located in the active site of the enzyme, and most of the residues involved in the reaction with these two ligands are positioned at the catalytic domain II, which in turn can have a ruinous impact on the enzymatic function of the protein. It deserves to mention that 2-Acetoxy-6-p-hydroxybenzoyl jaeschkeanadiol made two hydrophobic interactions with the catalytic dyad comprising Cys145 and His41, which shows this natural compound possesses a high potential in disturbing the deprotonating of Cys145 by His41, leading to a disruption of nucleophilic attack of Cys145 to polyproteins and their cleavage, and finally, disturbing the enzymatic function of the viral receptor. The other two natural compounds Culifolin and Teferin have bound to the active site with the binding energies of -8.2 and -7.9 kcal/mol, correspondingly, and with the amino acid residues of in the catalytic domain II (Figure 3 E&F).

These interactions can highlight their potent antiviral properties against COVID-19 M^pro.

Molecular dynamics simulations results

Molecular dynamics simulation is a highly effective approach for investigating the movements and interactions of macromolecules at both the molecular and atomic scales [7]. To investigate the effect of protein-ligands interactions on the protein structure, all systems consisting of the M^pro and ligands (with the lowest binding energy derived from docking results) were subjected to 150 nanoseconds of molecular dynamics simulation. The RMSD model can be used to identify when the system has reached equilibrium. The RMSD results can also be crucial in determining the appropriate simulation time for a target molecule [7]. As can be seen in the Figure 2, it can be asserted that 150 ns simulation period is sufficient for such a system because the system has gained its equilibration state during this period. Concerning RMSD analyses in complex with ligands, the complexes composed of (+)-Pinoresinol, 6 Cinnamylchrysin, 2-Acetoxy-6-p-hydroxybenzoyl jaeschkeanadiol, and 5,4'-Dihydroxy-3,3'-dimethoxy-2-prenyl-E stilbene have shown the highest values in terms of RMSD compared to the free protein. Furthermore, RMSD values for systems containing 2-Acetoxy-6-p-hydroxybenzoyl jaeschkeanadiol, (+)-Pinoresinol, and 6-Cinnamylchrysin ligands have indicated the highest instability in the enzyme building block (Figure 4).

Root Mean Square Fluctuation (RMSF) is an analytical method to evaluate the residues fluctuations in a protein. Essentially, RMSF is the standard deviation of each residue relative to the protein total displacement. By assessing the instability of amino acids, RMSF values will be higher. This value may vary for each amino acid depending on the degree of ligand-enzyme interaction. The attachment of all ligands to the M^pro enzyme led to changes the RMSF values of protein amino acids in different sections (Figure 5). In the protein-(+)-Pinoresinol complex, the RMSF values in amino acids 8-12, 24, 47, 72, and 83 (in domain I), and 97-101, 127, 154, and 177-178 (in domain II) have increased. However, the RMSF values in the location of amino acids 19, 35-33, 42-37, 45, and 50 (in domain I and in the active site), 60, 64, 69, 76, 80, and 93 (in domain I), 100, 102, 123, 137-141, and 143-146 (in domain II in the active site), 170-168, 190-186 (in the loop connecting domain II to III), 203-200, 217-214, 225-222, 231-227, 240-236, 243-241, 262-258, 266-264, 273-270, and 292-284 have decreased. This reduction may indicate the destructive effect of this compound on domain II, the active site, and also on the C-terminal of the protease enzyme (Figure 5A). In the complex of 6-Cinnamylchrysin with the protease, the RMSF values in amino acids 10-6, 12, 24, 43, 52-48, 64, 76, and 83 (in domain I), 97, 107-105, 117-113, 134-125, 139-137, 147-144, and 169-168 (in domain I), 184-179, 191-188 (in the loop connecting domain II to III), 198-194, 202-200, 232-225, 254-251, 256, 277, 290, and 294 (in domain III) have increased. However, these values have decreased in the location of residues 33, 118, 142, 155, 167-164, 217-214, 222, 239-236, 244, 275-273, and 288-282 (Figure 5B). Based on the RMSF data, it appears that the ligands attachment to the main protease enzyme of SARS-CoV-2 has significantly altered the RMSF values of the amino acids that play a significant role in the activity of this enzyme. This event can practically disrupt the protease enzyme's activity and prevent virus replication.

Regarding the protein in complex with 2-Acetoxy-6-p-hydroxybenzoyl jaeschkeanadiol, it can be mentioned that the RMSF values have mainly experienced fluctuations in the amino acids present in I and II domains (particularly the catalytic motif). These values, along with some of the amino acids present in the C-terminal (such as 255-248, 257, 270, and 298-291), have shown an increase in the amino acids with numbers 15-13, 37, 47, 47-72 (located in domain I), 107-105, 113-110, 135-125, 155-153 (located in region II), as well as 186-181 and 192-188 (present in the connecting loop between domains II and III). In addition, some decreases in the RMSF values of the amino acids present in the catalytic regions, such as 40-43 (in domain I) and 149-144 (in domain II), along with significant decreases in the amino acids located in the C-terminal, can contribute to the significant effect of this compound on the catalytic site, leading to disruption in enzyme function (Figure 5C). Furthermore, the effect on the C-terminal of the enzyme can affect its dimerization and the formation of its active form [7]. As for the RMSF plot of the enzyme in complex with 5,4'-Dihydroxy-3,3'-dimethoxy-2-prenyl-E-stilbene, it can be asserted that the binding of this natural compound to the enzymatic target leads to a decrease in the RMSF values of some amino acids in domain I, as well as some others present in the catalytic domain II, and the connecting loop between region II and III. In addition, it causes an increase in these values in some regions present in the active site of around domains I and II (Figure 5D).

   As for the RMSF of the other two remaining complexes comprising Culifolin and Teferin, the RMSF plots reveal that the most of the fluctuations in the RMSF magnitudes have occurred around amino acid residues present in the catalytic domains I and II, in particular in the latter system (Figures 5E&F). In general, the RMSF data suggest that the binding of ligands to the main protease of COVID-19, significantly alters the RMSF values of the amino acids that play a significant role in the catalytic function of this enzyme, which it can practically disrupt the enzyme’s activity and prevent viral replication. The Rg (radius of gyration) stands for of the total average dimension of a protein and is known as a protein stability indicator. The radius of gyration can also indicate the compactness of the protein. Any increase or decrease in this parameter can indicate a collapse or compression of the structure of the protein of interest. The compactness of enzyme structures can have a significant impact on the quality of substrate binding to the active site [7]. Compared to the free protein, the Rg values associated with the protein in complex with ligands decreased significantly, which could indicate compression of the enzyme structure after attachment to the ligands. However, it should be noted that the decrease in Rg value in the complex containing the Teferin was slightly less. Regarding the fluctuations in Rg value attributed to the alternating compactness in the enzyme structure, it could be noted that the highest fluctuations in Rg values ​​are related to complexes composed of (+)-Pinoresinol, 6-Cinnamylchrysin, 2-Acetoxy-6-p-hydroxybenzoyl jaeschkeanadiol, and 5,4'-Dihydroxy-3,3'-dimethoxy-2-prenyl-E-stilbene (Figure 6). Consequently, these notable fluctuations observed in the enzyme can provoke structural instability and disturb its enzymatic function.

   The hypothetical surface of solvent surrounding a protein, which comes into contact with the protein through Van der Waals forces, is known as the solvent-accessible surface area (SASA) of the protein [7]. SASA analysis was performed for protein-ligand complexes under investigation. As shown in Figure 7, the SASA values for all complexes have experienced a decreasing trend. This event could be due to the compression of the enzyme structure or the closure of water entry ports in the protein cavities, preventing water penetration into the interior of the enzyme. These results agreed with the Rg findings, and both emphasized the compression of the protein structure in an aqueous environment.

   Principal component analysis (PCA) is a method that can be employed to identify the primary constituents of protein movement. The essential movements that proteins undergo are critical for their structure and function. Any alterations in these movements can cause disturbances in the protein's activity [7]. The 2D patterns of M^pro movements were obtained separately and under conditions of ligand binding by PCA analysis, and these findings have been illustrated in Figure 8. According to Figure 8A, the free protein contains motion clusters with a moderately large motion range between -9 and 11 nanometers. It can be noted that almost all natural compounds under study have changed both the main protein movement pattern and its motion range, which could indicate the ability of these compounds to create instability and disruption in enzyme structure and function. However, the most significant shape change in the paradigm and the limitation of enzyme movements are related to complexes containing (+)-Pinoresinol, 2-Acetoxy-6-p-hydroxybenzoyl jaeschkeanadiol, 6-Cinnamylchrysin, and 5,4'-Dihydroxy-3,3'-dimethoxy-2-prenyl-E-stilbene, respectively. These limitations in the protein motion range (consistent with Rg findings) show that the three-dimensional viral enzyme structure is contracted or wrinkled in the vicinity of natural inhibitors. Such compressions in the structure of enzyme can hinder its ability to interact accurately with the substrate.

Pharmacokinetic properties and toxicity

The pharmacokinetic properties of six selected natural compounds were finally evaluated using the SwissADME server (Tables 2 and 3).

As the data shows, all of the natural compounds had excellent and desirable biocompatibility. Moreover, three compounds, (+)-pinoresinol and 6-Cinnamylchrysin, and 2-Acetoxy-6-p-hydroxybenzoyl jaeschkeanadiol, exhibited better pharmacokinetic properties (such as drug-likeness and pharmacological properties) compared to the other compounds. Furthermore, all of natural compounds except for 5, 4'-Dihydroxy-3,3'-dimethoxy-2-prenyl-E-stilbene possessed low or moderate toxicity. Based on this information, it can be mentioned that (+)-pinoresinol, 6-Cinnamylchrysin and 2-Acetoxy-6-p-hydroxybenzoyl jaeschkeanadiol did not show any drug-drug interactions with other drugs. However, due to the inhibitory impact of 6-Cinnamylchrysin on CYP2C9, it can be noted that this compound may have common interactions with the metabolism of some common non-steroidal anti-inflammatory drugs and sulfonylureas (Cytochrome P450 2C9-CYP2C9).

Conclusion

This study utilized molecular docking and molecular dynamics simulations to examine the potential of natural compounds found in propolis to act as antiviral agents against the COVID-19 main protease, 3CL^pro.

Results from the docking analysis identified six natural compounds, named including 6 Cinnamylchrysin, (+)-Pinoresinol, Culifolin, 5,4'-Dihydroxy-3,3'-dimethoxy-2-prenyl-E-stilbene, Teferin, and 2-Acetoxy-6-p-hydroxybenzoyl jaeschkeanadiol, with the strongest binding energies to the COVID-19 M^pro.

Further analysis was conducted on the pharmacokinetic properties of these compounds, and molecular dynamics simulations were performed to study the fluctuations in the protein's structure.

The simulations revealed that two of these compounds, (+)-Pinoresinol and 6 Cinnamylchrysin, caused the highest fluctuations in the protein's structure, which can pose negative effects on the enzymatic function of the viral target. However, these two compounds were found to have good pharmacokinetic properties and low toxicity, and they were able to cross the cell membrane, which is essential for them to access the viral protease. The study suggests that further research is needed to confirm the effectiveness of (+)-Pinoresinol and 6-Cinnamylchrysin as potential inhibitors of COVID-19 M^pro.

Disclosure Statement

No potential conflict of interest was reported by the authors.

Funding

The authors received no specific funding for this research work.

Authors' contributions

All authors contributed to data analysis, drafting, and revising of the paper and agreed to be responsible for all the aspects of this work.

Conflict of interest

All authors declare that there is no conflict of interest in this manuscript.

Acknowledgements

The authors would like to thank Imam Ali University for its support.

Orcid:

Armin Zarei: https://orcid.org/0000-0003-0761-1085

Rasool Amirkhani: https://orcid.org/0000-0002-2330-4598

Mahdi Gholampour: https://orcid.org/0000-0002-9553-4894

Hassan Tavakoli: https://orcid.org/0000-0002-9037-8300

Ali Ramazani: https://orcid.org/0000-0003-3072-7924

--------------------------------------------------------------------------

--------------------------------------------------------------------------

Copyright © 2023 by SPC (Sami Publishing Company) + is an open access article distributed under the Creative Commons Attribution License(CC BY)  license  (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.