In silico investigation of potential inhibitors to main protease and spike protein of SARS-CoV-2 in propolis.

Docking analysis of propolis's natural compound was successfully performed against SARS-CoV-2 main protease (Mpro) and spike protein subunit 2 (S2). Initially, the propolis's protein was screened using chromatography analysis and successfully identified 22 compounds in the propolis. Four compounds were further investigated, i.e., neoblavaisoflavone, methylophiopogonone A, 3'-Methoxydaidzin, and genistin. The binding affinity of 3'-Methoxydaidzin was -7.7 kcal/mol, which is similar to nelfinavir (control), while the others were -7.6 kcal/mol. However, we found the key residue of Glu A:166 in the methylophiopogonone A and genistin, even though the predicted binding energy slightly higher than nelfinavir. In contrast, the predicted binding affinity of neoblavaisoflavone, methylophiopogonone A, 3'-Methoxydaidzin, and genistin against S2 were -8.1, -8.2, -8.3, and -8.3 kcal/mol, respectively, which is far below of the control (pravastatin, -7.3 kcal/mol). Instead of conventional hydrogen bonding, the π bonding influenced the binding affinity against S2. The results reveal that this is the first report about methylophiopogonone A, 3'-Methoxydaidzin, and genistin as candidates for anti-viral agents. Those compounds can then be further explored and used as a parent backbone molecule to develop a new supplementation for preventing SARS-CoV-2 infections during COVID-19 outbreaks.

Introduction

Earlier this year, a novel coronavirus (SARS-CoV 2) triggered a coronavirus disease 19 (COVID-19) outbreak in Wuhan, China, and then spreading globally. The outbreak was announced as a Public Health Emergency of International Concern in January 2020 by World Health Organization (WHO). As of May 12, 2020, there are more than 4,324,739 confirmed cases of COVID-19, and at least 297,197 deaths reported, which indicates a severe public health threat [1]. Unlike MERS-CoV but similar to SARS-CoV [2,3], SARS-CoV 2 can cause human-to-human transmission, and its intermediate host leading to the current human infection is still under investigation [4,5].

Any concrete solutions in prevention, treatment, and post-infection are still lacking. So far, several drugs have been studied as COVID-19 treatments, such as nevilnafir [6], chloroquine [7], hydrochloroquine [8] and remdesivir [9]. Several researchers identified potential anti-viral agents to prevent SARS-CoV-2 infection [10,11]. Several FDA-approved drugs were evaluated as potential anti-viral agents. They can bind to SARS-CoV-2 RdRp and blocked its active site [10]. Remdesivir and chloroquine also have an inhibitor activity to SARS-CoV-2 infection [11]. Although these drug candidates have been massively studied, its currently undergoing the clinical evaluation phase and caution must be taken to avoid any negative effects [12]. Unfortunately, there is no apparent certainty which drug is the most effective in curing the SARS-CoV-2 infection.

It has been computationally determined that the SARS-CoV-2 has a similar mechanism and receptor to SARS, angiotensin-converting enzyme 2 (ACE2) [3]. SARS-CoV-2 uses its spike glycoprotein, which has similarities with the SARS-CoV spike protein structurally, to mediate membrane fusion. The conservation is 73% with most of the variability being in the protein's host cell interaction region. SARS-CoV-2 spike protein is a class I viral fusion protein, which required a protease cleavage to activate the membrane fusion. It has two subunit monomers, S1 for mediating attachment and S2 for membrane fusion, with a similar length of around 180 kDa. Any of them has been the main target for the neutralization antibody [13].

The second in silico docking model is the 3CLPRO main protease, which is responsible for controlling several major functions of the virus and has a highly conserved catalytic domain from the SARS virus [13]. The SARS-CoV-2 main protease (Mpro) is a key CoV enzyme that plays an essential role in mediating viral replication and transcription. When viruses enter the cell, RNA viruses take over the host cell to express structural and non-structural proteins, and Mpro cleaves the large viral polyprotein into a mature viral protein [14,15]. Some of its functions include the virus's replication processes, which makes it an ideal target for drug development [6]. Therefore, Mpro is also identified as an attractive target for inhibitory viral replication [16].

Propolis is commonly used as daily consumption or supplement due to its natural compounds with good pharmacological and pharmaceutical properties [17]. It is well-known propolis has been used as a folk medicine for 300 years B.C [18]. The propolis global world production was around 100 and 200 tons a year, and China produces propolis the most than other countries [19].

The consumption of propolis continues to this day, therefore, determining the bioavailability and safety of its compounds are essential. Furthermore, propolis has been known as an immune agent, works by increasing macrophages' performance [20,21]. Most of the propolis content belongs to flavonoid and polyphenol groups [22], wherein its anti-viral potential has been the subject of several reports [[23], [24], [25], [26], [27]] [[23], [24], [25], [26], [27]] [[23], [24], [25], [26], [27]]. Propolis is reported to inhibit varicella-zoster virus infection, herpesviruses, and HIV by inhibiting viral entry [[25], [26], [27]] [[25], [26], [27]] [[25], [26], [27]].

Several reports are investigating how flavonoids inhibit the infection of SARS-CoV-2 with its host. Some of the flavonoid from propolis such as kaempferol, quercetin, rutin, and galangin has been reported to inhibit SARS-CoV-2's Mpro [11]. Besides, rutin and kaempferol are reported to possess a binding affinity more than nevilnavir against SARS-CoV 2 Mpro [28,29]. Quercetin and rutin are also reported as valuable drug candidates for medicating SARS-CoV-2 based on a molecular docking result against SARS-CoV-2 Mpro and human receptor ACE-2 [30]. Both SARS-CoV and SARS-CoV-2 use ACE2 receptors to enter their host cell via their protein spike, which has a 76% homology [13]. This attachment only occurs in ACE2 but does not involve ACE [31]. Some reports indicate an increase in the high expression of ACE2 in patients with COVID-19 in the oral cavity and lung mucosa [[32], [33], [34]] [[32], [33], [34]] [[32], [33], [34]], and the blocking of spike protein (S2) is imperative to prevent the virus from penetrating the cell membrane. Another strategy targets the main protease of SARS-CoV-2 (Mpro) to prevent structural protein cleavage during the viral formation in the host cell.

This study investigated active propolis compounds' molecular interaction against SARS-CoV-2 main protease and spike protein by molecular docking approach. The molecular interactions between targets and inhibitors were analyzed thoroughly, and the mechanism for the difference in binding ability was also investigated. Hopefully, this information can be used as a reference in building synthetic compounds with better affinity or developing drug and vaccine candidates.

Materials and methods

Propolis screening

The propolis was supplied from PT Nano Herbaltama Internasional, South Tangerang, Indonesia. The preparation was done by adding 1 ml of propolis to 9 mL methanol in a 10 ml volumetric flask and then filtered with Minisart SRP15 Syringe Filter PTFE 0.2 μm (Sartorius, Germany). The filtered solution was added into a 2 mL vial for LC-MS/MS analysis. LC-MS/MS XEVO G2-XS Quadrupole-time-of-flight (The Waters, New Zealand) was used and separation was carried out using Acquity UPLC HSS-T3 column (100 × 2.1 mm, 1.7 μm). The column was maintained at 40 °C throughout the analysis, and the sample temperature was kept at 25 °C. The mobile phases used were 0.1% formic acid in water (A) and 0.1% formic acid in acetonitrile (B) at a flow rate of 0.6 ml/min. The elution gradient was set as follows: Starts at 1%, at 0.5 min 1% A, at 16 min 35% A, at 18 min 100% and 20 min 1%A. A total of 2 μl of the test solution was injected to screen flavonoids, and the chromatograph was recorded for 20 min. The mass scan range was set at 50–1200 Da with electrospray ionization (ESI) positive ion mode. The source temperature was 120 °C and the desolvation temperature was 550 °C. The capillary voltage was set as 2.0 kV and cone voltage was 40 V. The cone gas flow was set at 50 L/h and desolvation gas flow was 1000 L/h. The low collision energy was 6 eV and the high collision energy was 15–40 eV (ramp). 1 ng/mL solution of Leucine Enkephalin was used as external references (lock spray) for mass correction during acquisition. The Lock-Spray was infused at a flow rate of 10 μl/min and scan time was 0.1 s with an interval of 30 s.

Ligand retrieval and preparation

As much as 24 ligands were retrieved in 3D conformer with SDF formats from PubChem database (NCBI) (http://pubchem.ncbi.nlm.nih.gov). Of them, 22 ligands were the screened bioactive compound of propolis, and the other 2 were nelfinavir (C32H45N3O4S) and pravastatin (C23H36O7) with PubChem ID 64143 and 54687, respectively. Each ligand was optimized in molecular geometry, torsional barriers, and intermolecular-interaction geometry using MMFF94 partial forcefield in CHARMM and converted to a.pdb file using the Discovery Studio 2020 (BIOVIA)[35]. The ligand preparation only uses the propolis component's original ligand structure without considering the appropriate tautomeric forms of ligands, the specification of the chirality, and pH range addition.

Protein retrieval and preparation

Two proteins were used as the target for this study, SARS-CoV-2 main protease complex and SARS-CoV-2 S2 subunit protein. The crystal structures for those proteins were retrieved from Protein Data Bank (PDB) (www.rcsb.org) with the following PDB ID: 6LU7 for the complex of SARS-CoV-2 main protease (2.16 Å) and 6LXT for SARS-CoV-2 S2 subunit protein (2.90 Å). Each retrieved crystal structure was prepared using PyMOL version 2.0.4 program (The PyMOL Molecular Graphics System)[36] by removing the water molecule and ligand molecule. Prepared pdb file was then assessed for its Ramachandaran plot using RAMPAGE server (http://mordred.bioc.cam.ac.uk/~rapper/rampage.php) to check its readiness for docking analysis[37]. Binding site residues of each ligand were identified using the Computed Atlas of Surface Topography of Protein (CASTp)[38]. The binding pocket was shown as a red pocket of Mpro (Fig. 1(a)) and S2 (Fig. 1 (b)), whereas the key residue was presented in Table 1.

Fig. 1

The binding pocket of (a) SARS-CoV-2 Main Protease (Mpro) and (b) SARS-CoV-2 Spike Protein S2 unit.

Table 1

Grid Selection and targeted key residue.

Protein Target	Key residue	Center	Dimension (Å)
6LU7	Thr24, Thr25, Thr26, Leu27, His41, Thr45, Ser46, Met49, F139, Leu140, Asn141, Gly142, Ser143, C144, His163, Met165, Glu166, His172	X: 17.6528Y: 23.4131Z: 70.9264	X: 38.0184Y: 34.0608Z: 38.7272
6LXT	Leu922, Gln926, Asn928, Gly932, Lys933, Gln935, Asp936, Ser939, Ala942, Ser943, Gly946, Lys947, Gln949, Asp950, Val951, Asn953, Gln954, Asn955, Gln957, Ala958, Asn960, Thr961, Val963, Lys964, Gln965, Ser967, Ser968, Asn969, Phe970, Gly971, Ala972, Ser974, Asp1165, Leu1166, Gly1167, Asp1168, Ile1169, Ser1170, Gly1171, Ile1172, Asn1173, Ala1174, Ser1175, Val1176, Val1177, Asn1178, Ile1179, Gln1180, Lys1181, Glu1182, Ile1183, Asp1184, Arg1185, Asn1187, Glu1188, Val1189, Ala1190, Lys1191, Asn1192, Asn1194, Glu1195, Leu1197, Ile1198, Asp1199, Leu1200, Gln1201	X: −11.7071Y: 3.4251Z: −10.4247	X: 25.0000Y: 25.0000Z: 99.2687

Docking analysis, visualization, and validation

The prepared protein and ligands were loaded to PyRx (https://pyrx.sourceforge.io/home) and then prepared to be run in the built-in AutoDock VINA in the PyRx program[39].

Site-specific docking was performed with the grid selection parameter, as shown in Table 1, and the rest of the setting was left at default. The binding result table and the best model for each protein-ligand interaction were saved to be visualized. Docking visualization was done using Discovery Studio 2020 (BIOVIA) for its 3D conformation and 2D binding diagram [35]. Docking analysis and visualization were done in Microsoft Windows PC with Intel® Core™ i7-7700H CPU @3.6 GHz, 16 GB RAM with NVIDIA Titan XP GPU.

The docking analysis validation was inspired by the Shivanika et al. (2020) study [40], which used decoy ligand as a validator of the docking parameters and accompanied with positive control. Positive controls were selected based on the published research that confirms the chemical ligand binding (Table 2). The decoy substances were obtained from the DUD-E online server (http://dude.docking.org/) [41]. The docking parameters and process were identical to the previously mentioned process.

Table 2

Ligands for positive control.

No	Ligands name	Pubchem CID	Reference
Mpro (6LU7)
1	Caulerpin	5326018	[42]
2	Dexamethasone	5743	[43]
3	Dextromethorphan	5360696	[44]
4	Niclosamide	4477	[43]
5	Nitazoxanide	41684	[43]
6	Prednisolone	5755	[44]
S2 Protein (6LXT)
1	TGG	73178	[45]
2	Arbidol	131411	[45]
3	Luteolin	5280445	[45]
4	Quercetin	5280343	[45]

ADME and toxicity analysis

Absorption, Distribution, Metabolism, and Excretion (ADME) and Toxic analysis are important in supplement or drug design. Some properties, including molecular weight (MW), rule of five (RO5), acceptable and donatable hydrogen bonds, topological polar surface area (TPSA), polarizability, and estimated solubility (ESOL) were analyzed for compositions selected after docking. We consider the following properties for intestinal absorption, solubility levels, and Blood-Brain Barrier (BBB) penetration levels. The ADME parameters and pharmacokinetic properties' prediction was done by the SwissADME webserver (http://www.swissadme.ch/)[46]. While the toxicity information for acute effects was derived from PubChem, and some of the compounds that have no information in PubChem predicted using the consensus method in Toxicity Estimation Software Tool (T.E.S.T) was done [47].

Results and discussion

In this study, we successfully identified the compound of commercial propolis, which Indonesian commonly consumes. The profile of propolis compounds is shown in Table 3. A total of 22 compounds of polyphenol have been found in propolis, classified as flavonoid, phenolic acid, and phenylpropanoid. The flavonoid structure was built by fifteen carbon skeletons forming two benzene rings linked with pyrane heterocyclic ring. Moreover, three subgroups of flavonoids are known, such as flavones, flavonones, and flavonol, based on their position of hydroxyl bonds, carbonyl bonds, and hydrogen bonds [48]. 15 compounds identified as flavonoid in propolis at retention time 7.15, 7.16, 10.86, 11.2, 12.48, 12.51, 13.5, 13.55, 14.12, 14.34, 14.77, 15.08, 16.45 min were 3′-Methoxydaidzein, 3′-Methoxydaidzin, genistin, xanthomicrol, 3′,5,6,7-Tetrahydroxy-4′-methoxyisoflavone, methylophiopogonone A, 3′-Deoxysappanol, 3′4′,7-Trihydroxyflavanone, moslosooflavone, luteolin, 2′,6′-Dihydroxy-4′- methoxy-dihydrochalcone, chrysoeriol, jaceosidin, (3R)-7,2′,3′-Trihydroxy-4′- Methoxyisoflavanone and neobavaisoflavone, respectively. All flavonoid compounds found in propolis were classified into flavones group except for 3′-Deoxysappanol, which falls into the category flavonols. The rest were classified into phenolic acid, i.e., isoferulic acid, dimethylcaffeic acid, cinnamic acid, caffeic acid, benzyl caffeate phenylpropanoid (2,5- Dimethyl-7-hydroxychromone and isoaloeresin D). Isoferulic acid, cinnamic acid, caffeic acid, and dimethylcaffeic acid were classified into hydroxycinnamic acid sub-group, while benzyl caffeate was hydroxybenzoic acid [49].

Table 3

The identified compounds of propolis.

Compound	PubChem ID	Formula	Structure	Molecular ion (m/z)
Cinnamic acid	444539	C9H8O2	Image 2	149.0597
Caffeic acid	689043	C9H8O4	Image 3	181.0498
2,5-Dimethyl-7-hydroxychromone	5316891	C11H10O3	Image 4	191.0710
Isoferulic acid	736186	C10H10O4	Image 5	195.0648
Dimethylcaffeic acid	717531	C11H12O4	Image 6	209.0816
Benzyl caffeate	5919576	C16H14O4	Image 7	271.0960
3′4′,7-Trihydroxyflavanone	3496769	C15H12O5	Image 8	273.0768
2′,6′-Dihydroxy-4′ methoxydihydrochalcone	169676	C16H16O4	Image 9	273.1150
3′-Methoxydaidzein	5319422	C16H12O5	Image 10	285.0751
Luteolin	5280445	C15H10O6	Image 11	287.0553
3′-Deoxysappanol	13846660	C16H16O5	Image 12	289.1081
Moslosooflavone	188316	C17H14O5	Image 13	299.0922
Chrysoeriol	5280666	C16H12O6	Image 14	301.0705
(3R)-7,2′,3′-Trihydroxy-4′-methoxyisoflavanone	14353662	C16H14O6	Image 15	303.0862
3′,5,6,7-Tetrahydroxy-4′-methoxyisoflavone	10543410	C16H12O7	Image 16	317.0668
Neobavaisoflavone	5320053	C20H18O4	Image 17	323.1287
Jaceosidin	5379096	C17H14O7	Image 18	331.0820
Methylophiopogonone_A	10065830	C19H16O6	Image 19	341.1028
Xanthomicrol	73207	C18H16O7	Image 20	345.0969
3′-Methoxydaidzin	10527347	C22H22O10	Image 21	447.1278
Genistin	5281377	C21H20O10	Image 22	455.0953
Isoaloeresin D	76332505	C29H32O11	Image 23	595.1584

All the identified propolis compounds were docked against two proteins target: Mpro and S2. Several positive controls were also docked to the protein target. Additionally, a subset of docking decoys was added to validate the docking parameters further. Our docking procedure was validated using decoy ligands similar to nelfinavir which were obtained from DUD-E online server.

There were 51 decoy compounds matched with nelfinavir and docked against the active site of Mpro (6LU7) and S2 protein (6LXT). This was done to enhance ligand enrichment, which is essential to assess the docking procedure and eliminate false positives. The result's binding energies remained between −8.9 kcal/mol to −5.7 kcal/mol 33 decoy compounds out of 51 showed higher binding energy compared to the Mpro-Nelfinavir complex. In addition, 27 decoy molecules showed no binding with 6LU7. Hence, this confirms the docking efficiency and protocol. When conducting this study, no other studies discussed the potential interactions between ligands and the SARS-CoV-2 S2 protein (6LXT). Therefore, we are the first to take this approach.

Previously, we have tried docking against 6LU7 for several FDA-approved anti-viral drugs as its positive controls; Caulerpin, Dexamethasone, Dextromethorphan, Niclosamide, Nitazoxanide, and Prednisolone and displayed docking scores of −8.6, −6.9, −6.2, −7.4, −6.6, and −7.1 kcal/mol, respectively. We have docked several compounds (TGG, Arbidol, Luteolin, and Quercetin) based on Xiu et al. (2020) [45] against 6LXT as its positive controls. The binding affinity is −8.3, −6.7, −7.7, and −7.6 kcal/mol for TGG, Arbidol, Luteolin, and Quercetin, respectively. However, those compounds were not the type of FDA-approved anti-viral drugs. Pravastatin belongs to statins which are lipid-lowering, anti-inflammatory, and potentially antiretroviral drugs. Furthermore, on several cellular pathways, pravastatin can block the virus-host interaction system, which can act as an effective compound as opposed to H1N1 infection [50].

Nelfinavir is identified as a potential inhibitor against COVID-19 main protease, based on binding free energy calculations using molecular mechanics with generalized born and surface area solvation (MM/GBSA) model and solvated interaction energy methods [51]. Another study stated that nelfinavir has the second-highest docking score (−7.9 kcal/mol) among favipiravir, lopinavir, simvastatin, rosuvastatin, pravastatin, pitavastatin, lovastatin, fluvastatin, and atorvastatin (−5.8, −7.9, −7.0, −7.7, −6.6, −8.2, −7.4, −7.7, and −6.8 kcal/mol) against protease enzyme of SARS-CoV-2 (PDB: 6LU7). The binding sites involved are GLN189, GLU166, and CYS145, which were available in our docking result [52].

The binding affinity of each compound towards the protein target is shown in Table 4. The free energy of binding and the binding constant (Ki) are each other related to support the docking score data. So, the more negative the free energy of binding, the more favorable the binding, whereas the lower the value of Ki, the more favorable the binding is. They showed promising binding affinities towards Mpro and S2 with variable free binding energies ranging from −5.4 to −7.7 kcal/mol, and −5.3 to −8.3 kcal/mol, respectively. The wide range of binding affinity against the targeted protein may explain the existence of hydrophobic aromatic rings and hydrophilic hydroxyl groups of flavonoids [48].

Table 4

The predicted binding affinity energy of propolis compound against protein targets with its Ki value.

Compound	Sample name	Interaction with Mpro		Interaction with S2
		Binding affinity energy (kcal/mol)	Ki (M)	Binding affinity energy (kcal/mol)	Ki (M)
Cinnamic acid	P1	−5.4	0.9978	−5.3	0.9979
Caffeic acid	P2	−5.9	0.9976	−5.5	0.9978
2,5-Dimethyl-7-hydroxychromone	P3	−6.2	0.9975	−6.1	0.9975
Isoferulic acid	P4	−5.7	0.9977	−5.6	0.9977
Dimethylcaffeic acid	P5	−5.7	0.9977	−5.7	0.9977
Benzyl caffeate	P6	−7.1	0.9971	−6.5	0.9974
3′4′,7-Trihydroxyflavanone	P7	−7.5	0.9970	−7.6	0.9969
2′,6′-Dihydroxy-4′ methoxydihydrochalcone	P8	−6.9	0.9972	−6.7	0.9973
3′-Methoxydaidzein	P9	−7.3	0.9971	−7.6	0.9969
Luteolin	P10	−7.5	0.9970	−7.7	0.9969
3′-Deoxysappanol	P11	−7	0.9972	−7	0.9972
Moslosooflavone	P12	−7.4	0.9970	−7.4	0.9970
Chrysoeriol	P13	−7.3	0.9971	−7.7	0.9969
(3R)-7,2′,3′-Trihydroxy-4′-methoxyisoflavanone	P14	−7.1	0.9971	−7.5	0.9970
3′,5,6,7-Tetrahydroxy-4′-methoxyisoflavone	P15	−7.6	0.9969	−7.8	0.9969
Neobavaisoflavone	P16	−7.6	0.9969	−8.1	0.9967
Jaceosidin	P17	−7.2	0.9971	−7.3	0.9971
Methylophiopogonone_A	P18	−7.6	0.9969	−8.2	0.9967
Xanthomicrol	P19	−7.1	0.9971	−7	0.9972
3′-Methoxydaidzin	P20	−7.7	0.9969	−8.3	0.9967
Genistin	P21	−7.6	0.9969	−8.3	0.9967
Isoaloeresin D	P22	−7.4	0.9970	−7.8	0.9969
Nelfinavir	Control for Mpro	−7.7	0.9969	–
Pravastatin	Control for S2	–		−7.3	0.9971

There are no compounds found that have lower binding energy needed to bind Mpro compared to nelfinavir. The binding affinity of nelfinavir was −7.7 kcal/mol, while the other reports showed the binding energy of −7 kcal/mol [53,54]. We found the 3′-Methoxydaidzin (P20) has a similar binding affinity with nelfinavir against Mpro. In contrast, many compounds have lower binding energy compared to pravastatin (−7.3 kcal/mol) when docked against S2 except for cinnamic acid (P1), caffeic acid (P2), 2,5-Dimethyl-7- hydroxychromone (P3), isoferulic acid (P4), dimethylcaffeic acid (P5), benzyl caffeate (P6) and 2′,6′-Dihydroxy-4′ methoxydihydrochalcone (P8).

Interestingly, four compounds have binding affinity lower than −8 kcal/mol, which are neobavaisoflavone (P16, -8.1 kcal/mol), methylophiopogonone A (P18, -8.2 kcal/mol), 3′-Methoxydaidzin (P20, -8.3 kcal/mol), and genistin (P21, -8.3 kcal/mol). Those four compounds could also be embedded into the same binding pocket of Mpro with a similar orientation (Fig. 2 (a)). While S2, except for P18 and P20, shared the same binding pocket (Fig. 2 (b)). Therefore, those compounds (P16, P18, P20, and P21) were used in the molecular docking simulation against S2 to visualize the binding interaction. For comparison purposes, the same compounds were examined the binding interaction with Mpro.

Fig. 2

Binding orientation of propolis compounds with Mpro and S2 show they could be embedded into the same binding pocket of both protein receptors. (a) Binding orientation within Mpro, nelfinavir = red; neobavaisoflavone (P16) = yellow; methylophiopogonone A (P18) = blue; 3′-Methoxydaidzin (P20) = green; and genistin (P21) = orange. (b) Binding orientation within S2, pravastatin = red; neobavaisoflavone (P16) = yellow; methylophiopogonone A (P18) = blue; 3′-Methoxydaidzin (P20) = green; and genistin (P21) = orange. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 3 (a) shows the 2D visualization of the interaction between nelfinavir against Mpro. Two primary conventional hydrogen bonds were found in the nelfinavir, i.e., Glu A:166 and Gln A:189, whereby the residue Gln A:189 shares another bond as a carbon hydrogen bond. Hydrogen bonding gives a stabilization effect; thus, it plays an important role in biomolecular structures [55]. While other key residues interact through hydrophobic bonds (π-alkyl and alkyl) and electrostatic bonds (van der Waals). Hydrophobic bonds occur due to any interaction between hydrophobic amino acids with a polar solvent, i.e., water [56]. The orientation of nelfinavir and its orientation are shown in Fig. 3 (b).

Fig. 3

(a) Detail bonding of molecular interaction between nelfinavir and Mpro in 2D visualization, showing the presence of conventional hydrogen, carbon-hydrogen, π-Alkyl, and van der Waals bonding. (b) The 3D interaction of nelfinavir at the binding pocket of Mpro.

Glu A:166 and Gln A:189 were also found in the P18 (Fig. 4 (b)), P20 (Fig. 4 (c)), and P21 (Fig. 4 (d)) through O – H interaction, however, the hydrogen bonding only found in Glu A:166, whereas Gln A:189 forms van der Waals bond with Mpro. Contrary, only one hydrogen bonding (Tyr A:54 and hydroxyl group) was found in P16, as shown in Fig. 4 (a), without the presence of Glu A:166 and Gln A:189. It implies that there is a certain electron delocalization or charge transfer [57].

Fig. 4

2D and 3D molecular interaction between (a) P16, (b) P18, (c) P20, and (d) P21 against Mpro.

P16 has unique interactions with amino acid residues. As shown in Fig. 4(a), the presence of π interactions can stabilize the position of the P16 in the binding site. Sulfur functional groups in CYS A: 145 (C–S–H) form π-Sulfur and π-alkyl with aromatic ring ligands. The aromatic ring has a cloud of π-electron that can interact with a lone pair electron of S–H cysteine to form π-sulfur [58]. Furthermore, the σ-bond of C–S in cysteine may also interact with the π-electron cloud to form π-alkyl. Consequently, two interactions occurred; π-sulfur and π-alkyl; at key residue CYS A:145. A pair of key residue methionine; MET A:165 and MET A:49; form π-alkyl interactions with the same aromatic ring of P16. The two π-alkyls presence helps in the intercalation of ligands at the binding site [59]. The π-π T shape interaction gives a significant effect to binding energy because there are electrostatic, dispersion, induction, and exchange-repulsion contributions, while a π-electron cloud in an aromatic ring interacts with other π key residue in T shaped position [60]. The interaction of P16-Mpro is indicated by the interaction of the π system in an aromatic ring.

Hydrogen bond becomes a priority in determining the key amino acids, therefore, Glu A:166 might play an important role in intermolecular interaction between compound and Mpro. Hence, although P18, P20, and P21 have higher binding affinity than nelfinavir, these compounds are also candidates to inhibit the SARS-CoV-2 Mpro due to the existence of Glu A:166. A report showed a decrease of binding affinity energy to ~2 kcal/mol when replacing Glu166 with an alanine residue [54], indicates the importance of Glu A:166. Moreover, the other key residues were also found to form hydrogen bonding with Mpro, such as Ser A:144 and Leu A:141 in P18, Arg A:188, Thr A:190, Thr A:26 in P20, Thr A:190, Arg A:188, Thr A:26, Gln A:192, Asn A:142 in P21 as shown in Fig. 4.

Fig. 5 (a) shows the interaction between pravastatin and S2 in 2D and its orientation (Fig. 5 (b)). It showed six hydrogen bonding of key residues, i.e., Asn F:1194, Lys B:1191, Ser F:939, Arg E:1185, and Ser F:940. Among them, only Ser F:939 were found to be conserved in P16 and P21. In general, the hydrogen bonding in the P16, P18, P20 and P21 are less than 5, as shown in Fig. 6 (a) – (d). However, interestingly, the binding affinity is lower than pravastatin, as shown in Table 4. We believe the different binding affinity between pravastatin and P16, P18, P20, and P21 is due to the presence of π bonding.

Fig. 5

(a) Detail bonding of molecular interaction between pravastatin and S2 in 2D visualization, showing the presence of conventional hydrogen, Alkyl, and van der Waals bonding. (b) The 3D interaction of pravastatin at the binding pocket of S2.

Fig. 6

2D and 3D molecular interaction between (a) P16, (b) P18, (c) P20, and (d) P21 against S2.

π bond is non-covalent chemical bonds. For instance, π-alkyl (Lys B:1191) and π-cation (Arg E:1185) were found in P16, in addition to the hydrogen bonding of Asp B:1184 (Fig. 6 (a)). The aromatic ring of P16 binds to arginine in its position and forms of π-cation between the nitrogen (of the positively charged guanidine) and the π electron cloud of its aromatic ring. Electrostatic interactions dominate these bonds. Arginine itself is an amino acid with a positively charged guanidine side chain [61]. A report has reported that the non-covalent interaction of guanidine (arginine) -benzene iso-surface showed that the region of π-cation attraction is due to NH polarity, and the average π-cation strength of arginine has been reported to be −2.9 ± 1.4 kcal/mol [62]. The π-alkyl interactions are dominated by the dispersion strength between π and –CH groups of the amino acids, but these bonds are weaker than π-cations. It has been reported that CH4-π has interaction energy of-1.3 kcal/mol [60]. The hydrogen carbon bond is also observed in Fig. 6 (a). The carbon-hydrogen bond (CH ⋯ O) is weaker than the conventional hydrogen bond because CH O's average distance is longer than the conventional hydrogen bond (N–H ⋯ O, N–H ⋯ N, O–H ⋯ O). However, it is known that the carbon-hydrogen C–H O bond has an important role in the molecular recognition process [61,62]. The π-donor H is the interaction between the hydrogen bond and the π aromatic ring ligand system and has a significant effect on stabilizing the ligands' position at the binding site. All P16 aromatic rings were observed to form π interactions with residual amino acids, where these interactions have an important contribution to the binding energy.

Asp B:1199, Asn F:928, Gln F:935 formed hydrogen bonding with S2, while Glu B:1195 and Leu F:1197 was π-anion and π-σ bonding, respectively, as shown in Fig. 6 (b). The π-anion is defined as the electrostatic interaction between anions positioned above the aromatic ring where the binding energy is dominated by the electrostatic anion-induced polarity contribution [63,64]. Therefore, the side-chain carboxylic acid functional groups of glutamic acid interact with the π-electron cloud on the ligands' aromatic ring. The same aromatic ring is attached to the σ system in the hydrophobic side chain of leucine to form π-σ. It is thought that the two π interactions stabilize the conformational position of the ligands, while the other ligand side chains bond to the hydrogen at the binding site.

P20 has quite a lot of hydrogen bonding among the other compounds, which are Leu F:1200, Ile F:1198, Asn F:928, and Glu B:1195, as shown in Fig. 6 (c). Furthermore, the presence of π-bindings was also observed, such as π-cation, π-anion, π-σ, and π-alkyl. The lower binding affinity of P20 more likely due to the double bonding of key residue, for instance, Glu B:1195 interacts in hydrogen and π-anion bonding, while Asn F:928 forms hydrogen and π-σ bonding. The π-anion has binding energy that is proportional to π-cation, besides π-σ and π-alkyl have strong dispersion energy [60,65]. It is predicted that the π-interactions present a significant contribution to the ligands' binding energy to their receptors. Therefore, the binding affinity of P20 is −8.3 kcal/mol, which is less than 1 kcal/mol compared to pravastatin.

It should be noted that P21 did not have any hydrogen bonding, as shown in Fig. 6 (d), however, the binding affinity of P21 is similar to P20, i.e., −8.3 kcal/mol. It believes the presence of amide-π stacked bonding (Asn B:1187) and π-alkyl bonding (Lys B:1191 and Ala F:1190) greatly influence the binding affinity of P21. An amide-π bond is formed between π surface of the amide bond, which overlaps with the π on the surface of the aromatic ring ligand [62]. This interaction is proportional to the hydrogen bond's strength based on the contribution of electrostatic, induction energy, and dispersion energy to the bond energy. It has been reported that almost all amide-π-accumulation interactions occur between the protein amide backbone group and the ligand aromatic ring [66]. The π-alkyls also contribute to energy binding through their dispersion forces. It has been reported that CH4-π has an interaction energy of −1.3 kcal/mol [60].

The distribution of π -electron of ligands can interact with positive and negative charges of the amino acid side chain, also hydrophobic side chains, the π system as well as amida bonds, therefore forming π -cation, π -anion, π -σ, π- π interactions, and π -amide respectively. This bonding occurred when the ligand is face-to-face or face to edge position with key residue at the binding site. Apart from hydrogen, pi-bonds are important in biological structure [59,62,67].

Several studies have been reported that the amide-π interaction is important in structural biology [[68], [69], [70], [71], [72], [73]] [[68], [69], [70], [71], [72], [73]] [[68], [69], [70], [71], [72], [73]]. Such as in acetylated lysine (KAc), which showed highly conserved aromatic residues in the bromodomain [72], while bromodomain is a protein domain in most acyltransferases and associated with chromatin structure regulation. Also, the amide-π interaction exists in human deoxyhemoglobin bind to its drug [74], also for the binding of other ligands to enzymes [75], furthermore, to be a selective binding for an organic solution model system.

The occurrence of π-bonding may cause its enthalpy/entropy on free binding energy be minimal. This interaction is unaffected by the solvation/desolvation process, which breaks water-mediated bonds and residue position rearrangement. π interaction has various binding confirmation which explains why there is various type of π interaction was detected between propolis' compound and the targeted protein (Mpro and S2), such as face-to-face (π-σ, π-cation, and π-π), parallel-displace and edge-to face (both occurring in π-π) [76].

Understanding hydrophobic interaction at molecular length scales may support a better understanding of protein folding, which remains unclear [77,78]. Hydrophobic amino acids are including Ala, Gly, Ile, Leu, Pro, Val, and Phe, which have hydrocarbon sidechains due to their non-polar chemistry [[79], [80], [81]] [[79], [80], [81]] [[79], [80], [81]]. In general, nelfinavir – Mpro interaction involves more hydrophilic amino acids (Glu, Gln, Cys, His and Asp) rather than hydrophobic amino acids (Leu and Val) also amphipathic amino acid (Met), which makes this bonding more hydrophilic than hydrophobic. The four potential inhibitors (P16, P18, P20, and P21) had a similar hydrophilic profile as nelfinavir - Mpro. Whereas pravastatin – S2 also interaction involves hydrophilic amino acids (Glu, Gln, Asp, Asn, Ser, and Arg) rather than hydrophobic amino acids (Ala) also amphipathic amino acid (Lys). This bonding is more hydrophilic rather than hydrophobic, also then nelfinavir – Mpro hydrophilic interaction. Among the four potential compounds (P16, P18, P20, and P21), P21 has the highest hydrophilic profile similarity with pravastatin – S2 interaction. These profiles crucially may affect herbal compound's ability to dissolve in the gastrointestinal fluids [82]. It is further explained in Table 5, which shows each compound's lipophilicity was according to XLOGP3. However, how the hydrogen bond at the identified ligand-protein interface relates to the hydrophobic interactions, which may increase the binding affinity has not been properly elucidated.

Table 5

ADME properties prediction of all propolis’ compounds.

Sample	MWa(g/mol)	RO5b	#H-bondacceptorsc	#H-bonddonorsd	TPSAe(Åb)	XLOGP3f	ESOL(mol/l)g	log Kp(cm/s)h	LD50(mg/kg)i
P16	322.35	0	4	2	70.67	4.40	1.24E-05	−5.14	768.89
P18	340.33	0	6	2	89.13	3.86	1.89E-05	−5.64	206
P20	446.40	0	10	5	159.05	0.64	6,98E-04	−7,05	934.65
P21	432.38	1	10	6	170.05	0.86	6.60E-04	−8.33	2000

Molecular Weight (acceptable range 130–725 gm/mol).

Rule of five Lipinski rules.

Acceptable H-bonds (acceptable range 2–20).

Donatable H-bonds (acceptable range 0–6).

Topological Polar Surface Area (TPSA < 140 Åb good intestinal absorption) (TPSA < 70 Åb good brain penetration).

Lipophilicity descriptor (Acceptable range -2 – 6.5).

Estimated solubility in water.

Skin permeant (Acceptable range (−8) - (−1)).

Predicted Oral Rat LD50 (Acute Effects).

Those potential compounds of the propolis (P16, P18, P20, and P21) were evaluated for their ADMET physicochemical properties prediction using SwissADME and toxicity using T.E.S.T, as shown in Table 5. P21 has the highest MW among the compounds, while P16 is the lowest. Indeed, the MW of all compounds within the acceptable range. All of these compounds obey Lipinski rule, except for P21, as indicated by the RO5 value higher than 0. For the H-bond acceptor and donor, all compounds showed the value in the acceptable range, even though H-bond donor P21 is in the edge, i.e., 6. Each compound's pharmacokinetics was also evaluated, which reveals P16 and P18 may be absorbed by gastrointestinal organ due to topological polar surface area (TPSA) is less than 140 Å2. The lipophilicity of compound was according to XLOGP3 value and showing its polarity. The lower the value, the higher the polarity. Only P21 was out of the acceptable range of skin permeant (log Kp) with a value of −8.33 cm/s. Toxicity analysis was done to predict each compound's LD50 based on Hodge and Terner (1995). P16, P20, and P21 belong to the slightly toxic drug category (500–5000 mg/kg) with LD50 of 768.89, 934.65 and 2000 mg/kg, respectively. Whereby P18 is in the moderately toxic drug category (50–500 mg/kg) with LD50 of 206 mg/kg.

To the best of the authors' knowledge, this study is the first to report the potential of 3′-Methoxydaidzin, genistin, methylophiopogonone A, and neobavaisoflavone to inhibit virus infection. 3′-Methoxydaidzin of propolis as the best binding ligand predicted for both protein targets; SARS-CoV-2 Mpro and S2. Blocking its spike protein may prevent SARS-CoV-2 bonding with the ACE2 receptor. Moreover, 3′-Methoxydaidzin can inhibit the extensive proteolytic activity of viral polyproteins [[83], [84], [85]] [[83], [84], [85]] [[83], [84], [85]]. 3′-Methoxydaidzin has no profound information except for its neuroprotective effect [86]. Interestingly, this chemical compound is extracted from Radix puerariae or kudzu plant commonly used as traditional Chinese medicine to treat diarrhea, acute dysentery, deafness, and cardiovascular diseases [87]. Genistin has evaluated to be anti-microbial, anti-cancer, anti-oxidative, cardioprotector, anti-apoptotic, neuroprotector, hepatoprotective [88], and anti-HIV [89]. Genistin is also found in kudzu plants, and a previous study showed that this compound could block the tyrosinase enzyme in its active site [[87], [88], [89], [90]] [[87], [88], [89], [90]] [[87], [88], [89], [90]].

On the other hand, methylophiopogonone A presents an anti-inflammatory activity [91,92]. Some of the publications about methylophiopogonone A are related to antioxidant activity in the form of herbal medicine supplementation [[87], [88], [89], [90]] [[87], [88], [89], [90]] [[87], [88], [89], [90]]. Lipinski [93] indicates a higher bioavailability is the compound's character related to lower molecular weight, lower H-bond capacity, and low lipophilicity. Furthermore, neoblavaisoflavone, methylophiopogonone A, 3′-Methoxydaidzin, and genistin exhibit the potential to be good for oral bioavailability. The median lethal dose (LD50) is a standard measurement for acute toxicity that predicts the adverse effect of substances in a short period after exposure [94]. No compounds of propolis in our study indicate highly toxic or extremely toxic. Propolis and its constituent are well tolerated nontoxic based on clinical and research reports [17,95,96].

Conclusion

In this study, we examined the propolis compounds to inhibit the SARS-CoV-2 Mpro and S2 using a molecular docking simulation. The propolis examined in this study has 22 compounds, which can be categorized as flavonoid, phenolic acid, and phenylpropanoid. The candidate to inhibit the SARS-CoV-2 Mpro was methylophiopogonone A, 3′-Methoxydaidzin, and genistin. Although the binding affinity energy of those compounds is not lower than nelfinavir (control), the existence of GLU A:166 in those compounds is essential to inhibit SARS-CoV-2 Mpro. In addition to GLU A:166, many other hydrogen bonds were observed in those compounds. For the blocking of SARS-CoV-2 S2, four propolis compounds (neoblavaisoflavone, methylophiopogonone A, 3′-Methoxydaidzin, and genistin) have binding affinity energy lower than the pravastatin (control), makes those compounds a candidate to inhibit SARS-CoV-2 S2. Although only neoblavaisoflavone and methylophiopogonone A has the ability to absorbed by gastrointestinal organ, the main purpose of propolis is used as a supplement instead of a drug. Moreover, the finding on methylophiopogonone A, 3′-Methoxydaidzin, and genistin as a candidate for anti-viral reported for first time in this study. Thus, those compounds can be further explored and could be used as a parent backbone molecule to develop new supplementation for preventing SARS-CoV-2 infections during COVID-19 outbreaks.

CRediT author statement

Azza Hanif Harisna: Conceptualization, Methodology, Software, Writing - Original Draft. Rizky Nurdiansyah: Conceptualization, Methodology, Software, Writing - Original Draft. Putri Hawa Syaifie: Formal analysis, Investigation, Writing - Original Draft, Project administration. Dwi Wahyu Nugroho: Validation, Formal analysis, Investigation, Writing - Original Draft. Kurniawan Eko Saputro: Formal analysis, Investigation, Data Curation, Writing - Original Draft. Firdayani: Formal analysis. Chandra Dwi Prakoso: Formal analysis. Nurul Taufiqu Rochman: Resources. Nurwenda Novan Maulana: Resources. Alfian Noviyanto: Visualization, Supervision, Writing - Review & Editing. Etik Mardliyati: Resources.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.