Dietary polyphenols mitigate SARS-CoV-2 main protease (Mpro)-Molecular dynamics, molecular mechanics, and density functional theory investigations.

The recent evolution of the SARS-like Coronavirus has ravaged the world. The deadly virus has claimed over millions of lives across the world and hence highlights the need to develop effective therapeutic drugs to contain the disease posed by this parasite. In this study, the inhibitory potential of fifty (50) dietary polyphenols against Coronavirus (SARS-CoV-2) main protease (Mpro) was conducted using the Autodock Vina Molecular docking tool. In the virtual screening process, the binding affinity of Remdesivir (-7.7 kcal/mol) currently used to treat COVID-19 patients was set as the cut-off value to screen out less probable inhibitors. Ellagic acid, Kievitone, and Punicalin were the only promising ligands with binding affinities (-8.9 kcal/mol, -8.0 kcal/mol and -7.9 kcal/mol respectively) lower than the set cut-off value. Furthermore, we validated Ellagic acid and Kievitone efficacy by subjecting them to molecular dynamics simulation and further stability was assessed at the molecular mechanics and quantum levels. The overall analysis indicates both compounds demonstrate higher stability and inhibitory potential to bind to the crucial His41 and Cys145 catalytic dyad of Mpro than the standard drug. However, further analysis of punicalin after evaluating its docking score was not conducted as the ligand pharmacokinetics properties suggests it could pose serious adverse effect to the health of participants in clinical trials. Hence, we employed a more safe approach by filtering out the compound during this study. Conclusively, while Ellagic acid and kievitone polyphenolic compounds have been demonstrated to be promising under this in silico research, further studies are needed to substantiate their clinical relevance.

Introduction

In this in silico-based research, we used ADMET profile, physicochemical properties analysis and molecular docking simulation to investigate some class of dietary polyphenols that could inhibit the viral growth of SAR-CoV-2 through the blocking of SARS-CoV-2 main protease 3CL (Mpro). We further subjected the selected compounds to 20 ns all atoms molecular dynamics simulation to investigate their stability at the Mpro active pocket. Molecular mechanics MMPBSA method was used to calculate the free energy of the residues involved in the interaction of the selected compounds with Mpro while Quantum mechanical-Density Functional Theory (HOMO and LUMO) was exploited to calculate the energy donating and accepting potential of these compounds. Our aim is to investigate the compounds that could emerge and developed as potent nutraceuticals for the management/treatment of SARS-CoV-2 mediated COVID-19.

Several reports of respiratory illness have been reported from Wuhan, Hubei province of China in December 2019, with confirmed evidence that the infections were caused by a novel coronavirus which was subsequently named SARS-COV-2 by World Health Organization [12,16,23,52]. The pathogenesis of SARS-COV-2 leads to a medical complication named COVID-19. This pathogen is a positive-sense single-stranded RNA virus from the genus β-coronavirus. The outbreak of this disease has become a pandemic to the world and a serious health challenge to the public health organization, owing to an historical record of about two million death cases. It was primarily believed to have been the first host to infect human, thorough investigation has established SARS-COV-2 to be closely related to SARS-CoV which occurs in 2003 and MERS-CoV in 2012 [3,5,17,38]. Out of the different therapeutic drugs and vaccines that have been proposed to reduce high increase or spread of this coronavirus disease, only the broad spectrum nucleotide analogue prodrug remdesivir (Velour) was the first and only approved drug by the FDA for the treatment of SAR-CoV-2-induced COVID-19 for hospitalized adults and children of age 12 years and above weighing at least 40 kg. In a research where 541 patients were assigned to Remdesivir and 521 to placebo, it was reported that Remdesivir was superior to placebo through the shortening of recovery time in COVID-19 infected hospitalized adults with evidence of lower respiratory tract infection [5], therefore, an urgent need to find effective natural anti‐SARS‐CoV‐2 measures for the treatment/management of COVID‐19 to alleviate the current global public health threat is the order of the day in the scientific world.

Persistent investigation of SARS-COV-2 genomics organization reveals that the viral genome has approximately 30 kb with at least 14 open reading frames which is released into the cytoplasm [21,30]. The ORF1a and ORF1ab encodes all non-structural proteins, the remnant genes encode for the structural proteins; spike (S) protein, envelope (E) protein, membrane (M) protein and nucleocapsid (N) protein, and nine accessory factors which together function in the virion formation. Mostly, the nonstructural proteins are important for the replication and transcriptions of the virus and it has been observed that a slippery sequence UUUAAAC is present at the junction between ORF1a and ORF1ab in all coronaviruses. Translation begins at the end of the slippery sequence through a RNA-mediated ribosome frame shift. The ORF1a encodes the papain-like (PLpro) and 3CL-protease (Mpro) which cleaves the replicates polypeptide to encode 16 non-structural proteins (NSPs) [15]. It has been reported that 3CL-protease proteolytic activity is greater than the Papain-like protease by cleaving larger number (11) of sites within the polyprotein, hence it is regarded as the major proteolytic enzyme of SARS-COVs [19,48]. The individual homo-dimeric protomers of SARS-COV2 Cysteine-like protease (Mpro) consist of three different domains: domain I consist of sequences between 8-101 amino acid residues, domain II is found at the region of 102–184 residues and sequence between 201-303 is made up of domain III [26]. The domain I and II which contains the catalytic dyad (His41 and Cys145) at the formed cleft [9,26,43] are made up of β-sheet while the dominant structural motifs of domain III is α-helix [26]. Crucially, Cys145 residue act as a nucleophile and has been identified as the major key player in the proteolytic activity of Mpro [26]. Hence, targeting this residue has been identified as a major therapeutic strategy in the prevention of polyprotein processing and Viral maturation exhibited by the Mpro [53]. Interestingly, we have recently reported the essence of targeting relevant amino acid residue on protein targets in computer-aided drug discovery [7]. In 2020, Jin and collaborators employed the use of Computer-aided drug design (CADD) and X-ray crystallographic technique to elucidate the structure of Mpro bound with a mechanism-based inhibitor (N3) [26]. This discovery as thus paves the way for researchers to investigate potential Mpro inhibitors from natural products/phytochemicals and clinically approved unrelated COVID-19 on Mpro (for drug repurposing) and other protein targets using in silico approach [[1],2,10,18,40].

Polyphenols majorly synthesized from phenylpropanoid pathways by plants are one of the phytochemicals which has been scientifically supported by nutritional and medicinal evidences as established functional foods or nutraceutical that makes the lines between foods and drugs to fade. They are polymeric unit of phenyl with more than a double phenyl rings with single or additional hydroxyl substituents. Food sources such as black tea, walnut, apples, blueberries, vegetables, almonds, spinach, tomatoes, parsley, leguminous plants, soy beans, cummoet, pomegranate, cherries, persimmons, rooibus tea, flex-seed etc., housed different subclasses of polyphenols in different proportions [4,8]. Long consumption of diet rich in polyphenol potentially prevent against infections and non-communicable diseases such as lungs damage, neurodegenerative diseases, osteoporosis, pancreatitis and type-2 diabetes. Several scientific examinations have established that various polyphenol have anti-inflammatory, antioxidant, antiviral activities and also have effect in numerous cellular processes. Ellagitannins possess antiviral activities in particular against HIV infections [39] and manifest inhibitory effects on Epstein-Barr virus, HSV-1 and HSV-2. Ellagic acid and Kievitone have been found to possess antioxidant, anti-cancer and antiviral activities [20,22,44]. Pre-existing drugs have been used by scientist in both ligand-based and structural-based method to control the spread of COVID-19, however some of these drugs have been thus far had inappreciable effects on SARS-COV-2 which could be due to the suggestion that the viral proteins may have actively acquire mutations by single nucleotide polymorphisms (SNPs) and thus escape for being targeted by antiviral drugs. However, the use of natural product such as polyphenols may leverage this insurgence and could be a promiscuous therapeutic target against SARS-COV-2 proteins [51].

Materials and methods

Preparation of target protein

The crystal structure of the SARS-COV-2 Main protease PDB ID: 6lu7 [26] was harvested from the RCSB Protein Data Bank (http://www.rcsb.org). We validated the protein by checking the presence of missing hydrogens, side-chain abnormalities, iMproper bonds etc. The mechanism-based inhibitor (N3) which was co-crystallized with Mpro was removed from the protein before executing our docking protocols. The protein structures were inserted in Autodock Vina and the pdbqt format for each protein was generated [49].

Determination of (6LU7) Mpro active sites

Binding pocket, ligands interactions and all amino acid residues in the active site of SARS-CoV-2 Mpro were determined with Computed Atlas for Surface Topography of Proteins (CASTp) (http://sts.bioe.uic.edu/castp/index.html?2011) and Biovia Discovery Studio (2019). Obtained data were compared and validated with previously reported experimental data for SARS-CoV-2 proteins.

Ligands selection and preparation for docking

We retrieve the canonical smiles for 3-hydroxyflavone, Apigenin, Baicelin, Caffeic acid, Catechin, Curcumin, Cyanidin, Diadezin, Delphinidin, Ellagic acid, Epiafzellichin, Epicatechin, Epigallocatechin, Ferulic acid, Gallic acid, Gallocatechin, Genisterin, Glycitein, Hesperetin, Hesperidin, Isorharmetin, Kampferol, Kievitone, Lariciresinol, Malvidin, Medioresinol, myricetin, Naringenin, Naringin, Narirutin, Neohesparindin, Nobiletin, p-Coumatic acid, Pelargiridin, Peonidin, Petinidin, Pinoresinol, Protocatechuic, Punicalin, Punicalagin, Quercetin, Resveratrol, Rutinoside, Secisolariciresinol, Sesamin, Sinapic acid, Syrigaresinol, Tangeretin, Vallic acid and Woogin acid and also the standard (Remdesivir) from PUBCHEM (http://pubchem.ncbi.nlm.nih.gov). Cactus Translator Server (https://cactus.nci.nih.gov/translate/) was employed to generate pdb format for each compound canonical smiles. The compounds were optimized and inserted into Autodock Vina to obtain their respective pdbqt formats [49].

Molecular docking

Further analyses were carried out using Autodock 4.2 for identifying specific interaction involved in binding of molecules to the target receptor. Subsequently, all the proteins were preprocessed and Gesteiger charges were added to it using AutoDock MGL tools 1.5.6. A receptor grid box was generated after the active sites had been defined. Grid box dimension of 60×60×60 for 3cl Mpro with spacing of 0.375 Å centering on residue Cys145 and His41 was defined. The protein-ligand interactions were visualized by Discovery Studio and Pymol [11].

ADMET and drug-likeness analysis

An efficient and widely used computational server-Admetsar2 (https://lmmd.ecust.edu.cn/admetsar2/) was entrusted to perform ligand-based studies of the selected phytochemicals to ascertain the compound that may have the potential of qualifying to be a druggable lead compounds with the aim to reduce the viral pathogenesis. Diverse parameters under numerous rules were used to screen the pharmacokinetics and drug-likeness of these compounds [45,54].

Molecular dynamics simulation

Prior to the 20 ns production run of the three complexes (including the standard), GROMACS (GROningen MAchine for Chemical Simulations) Molecular dynamics operating environment [50] was used to prepare the files for the simple dynamics normal control group “6LU7” and for the complexes which include 6LU7-Ela, 6LU7-Kie, and 6LU7-Rem. We used CHARMS-36 force field and TIP3P GROMACS recommended water model for the protein topology. CGENFF web server tool was used to prepare the “.str” file for the ligand topology after which the topology files of the complexes were appropriately updated using appropriate python codes to manually include ligands topology. We ionized and neutralized the systems and then solvated them using the simple point charge-216 explicit water model (spc216.gro). Energy minimization was run for 100 ps (picoseconds) using steepest descent algorithm to establish a stable conformation for the system after which we equilibrated the system using Verlet algorithm from 0-310 K at 100 ps, 2fs (femtoseconds) time step for NVT, while Berenson algorithm at 2fs time step for 100 ps was used for NPT equilibration. The system was treated using the “trjconv” module to centralize and compact the protein while we also made sure the atoms does not jump out of the PBC (Periodic Boundary Condition). RMSD (Root Mean Square Deviation) was calculated and reported as mean ± SD using the 2004 frames for each 20 ns production step run for each complex using “gmx analyze” GROMACS module after which the H-bonds were also calculated using “xmgrace” module.

Molecular mechanics (MM-PBSA)

The Molecular Mechanics Poisson–Boltzmann surface area (MM-PBSA) method [31], implemented in the g_mmpbsa module [32] of GROMACS, was used to compute the relative binding free energy of each complexes for the last 5 ns trajectories. Xmgrace was then used to plot the graph spectrum of the individual contributing energy of the residues for the three complexes.

The Binding energies consist of the van der Waal forces, electrostatic interactions and the solvation energy while this solvation energy comprises of both polar and non-polar energies. The equation is defined below as:Where BE = Binding Energy, VDW =Van Der Waal, PSE = Polar Solvation Energy, and SASA=Solvent Accessible Surface Area Energy

Density functional theory (quantum mechanics)

The top two compounds (hits) from the docking screening and the standard (Remdesivir) were subjected to quantum mechanical calculation using density functional theory. The Gaussian 09 W program [13] was used for the calculations by optimizing the compounds' geometries at DFT/B3LYP/6–31G(d'p’) levels. The frontiers orbital energies, the highest occupied molecular orbital (HOMO), the lowest occupied molecular orbital (LUMO), energy gap, and the molecular electrostatic potential were computed in order to understand the electron acceptor and electron donor properties of the compounds. These also provide information about the chemical reactivity and stability of the compounds.

Results and discussion

Molecular docking experiment

Right here for the first time, we report the anti-SARS-CoV-2 potential of two polyphenols (Ellagic acid and Kievitone). We study the inhibitory potential of 50 dietary polyphenols against the Mpro target (6LU7) of SARS-COV2 using Autodock Vina. The binding affinity result of the 50 dietary polyphenols has been depicted in supplementary data (Table 1 ). The computational thermodynamic scoring function of the widely used Autodock Vina program reveals three top-ranking compounds-Ellagic acid, Kievitone and Punicalin with significant binding affinity of −8.9 kcal/mol (4 conventional H-bond and 3 hydrophobic bonds), −8.0 kcal/mol (8 H-bond and 3 hydrophobic bonds) and −7.9 kcal/mol (6 H-bonds and 1 hydrophobic interaction) respectively (see Table 1 and Fig. 1 ). Other compounds were filtered out based on the binding affinity of the standard drug-Remdesivir (−7.7 kcal/mol) which was used as cut-off to screen out less potent compounds (Table 1). Since lowest binding energy ligand represents highest binding affinity compound [47], Ellagic acid could be considered as the best docked compound in this study.

Table 1

indicates the binding affinity, Inhibition constant, hydrogen bonds and hydrophobic bonds interaction of dietary polyphenols hits against SARS-COV-2 main protease (6LU7).

Compounds	Targeted Proteins 6LU7
	Binding Affinities	Hydrogen Bonds interactions	Hydrophobic Bonds Interactions
Ellagic acid	−8.9 kcal/mol	His41, Asn142, Gly143, Glu166	Met49, Cys145
Kievitone	−8.0 kcal/mol	Asp187, Cys143, asn142, Glu166, His172, his163, Ser144	His41, Met49, Met165
Punicalin	−7.9 kcal/mol	Thr26, Leu141, Gly143, His164, Met165, Asp187, Arg188,	Cys145
Remdesivir	−7.7 kcal/mol	Glu166, Ser144, gly143, Cys145, asn142, His163	His41, His164, Gln165

Fig. 1

3D structure of SARS-COV2 Mpro and A. Ellagic acid B. Kievitone C. Punicalin D. Remdesivir conventional hydrogen bonds (Green), carbon hydrogen bonds (cyan) and hydrophobic (pink) 2D interactions with SARS-COV-2 Main protease residues.

Moreover, our result of conformational sampling of Ellagic acid and Kievitone against the Mpro molecular target also indicates that the phytocompounds binds efficiently with domain I and II of the receptor active site including the His41 and Cys145 catalytic dyad (Fig. 1). Although Punicalin only interacts with Cys145, however, the ligand depicted similar mode of binding when docked with the receptor. Interestingly, similar mode of binding of these phytochemicals was also observed in the standard drug (Remdesivir) which was also found interacting with the crucial catalytic dyad of Mpro (Fig. 1). Several investigations which support our findings have been reported in the studies of Gurung and collaborators, Das and colleagues, Joshi and Co-authors and Islam et al. [10,18,24,27]). We believe the competitive binding of our hit compounds with His41 and Cys145 could be a potent therapeutic measure that could inhibit the SARS-COV2 Mpro catalytic activity.

ADMET and drug-likeness

Computational drug-likeness test is done to facilitate the creation of new drug molecules, the molecules of which must comply with the criteria of the following five Lipinski's rule: molecular weight must not be greater than 500, number of hydrogen bond donors must not be greater than 5, number of hydrogen bond acceptors must not be more than 10 and lipophilicity (expressed in LogP) must not be greater than 5 [34]. Potential inhibitors with more than one violation of these rules are considered as molecules with poor absorption and low permeation properties [6]. The result of this study shown in Table 2 reveals the potential of Ellagic acid and Kievitone as an excellent drug-like compound since both compounds do not violate any of Lipinski rules. In contrast, Punicalin do not comply with Lipinski drug-likeness with three (3) violations of the physiochemical properties of the rule of five. Hence, there is high propensity that Ellagic acid and Kievitone could behave as a better drug-like molecule than the standard which disobeyed three properties from the rule of five.

Table 2

Indicates the physicochemical properties of selected hit compounds and the standard used in this virtual screening.

Compounds	Lipinski's Rule of 5
	Mol. wt.	Hydrogen donor	Hydrogen Acceptor	Log P	No. of Violation
Ellagic acid	302.19	4	8	0.79	0
Kievitone,	356.37	4	6	2.24	0
Punicalin	782.53	13	22	−0.07	3
Remdesivir	602.58	4	12	3.24	2

In silico ADMET analysis is used to assess the pharmacokinetic and toxicological properties of hit compounds before experimental procedures [25]. ADMET evaluation of compounds is crucial in determining whether a potential lead candidate is capable of passing through the clinical trials stage. Hence, it is imperative to perform ADMET test early during the drug discovery process in order to filter out small molecule inhibitors with unacceptable and poor pharmacological profile. Our result of pharmacokinetic predictions of Ellagic acid, Kievitone, Punicalin and the standard drug is represented in Table 3 .

Table 3

Represents the ADMET properties of lead compounds selected from the bulk of dietary polyphenols used in this virtual screening. C-1= Ellagic Acid; C-2= Kievitone; C-3 = Punicalin C-4= Remdesivir.

Absorption&Distribution	C-1	C-2	C-3	C-4
BBB (+/-)	–	–	–	+
p-glycoprotein substrate	–	–	–	+
HIA+	+	+	+	+
Metabolism
CYP450 2C19Inhibitor	–	+	–	–
CYP450 1A2Inhibitor	–	+	–	–
CYP450 3A4Inhibitor	–	–	–	–
CYP450 2C9Inhibitor	–	+	–	–
CYP450 2D6Inhibitor	–	–	–	–
Excretion
Biodegradation	–	–	–	–
Toxicity
AMES Mutagenesis	–	–	+	–
Acute Oral Toxicity	II	III	III	III
hERG Inhibition	–	–	+	–
Carcinogenicity	–	–	–	–
Hepatotoxicity	+	+	+	+

Human intestinal absorption (HIA), substrate of P-glycoprotein and permeation through the blood-brain barrier (BBB) have been identified for the prediction of absorption and distribution properties of our drug candidates (Table 3). The positive value of HIA depicted by the selected compounds including the standard may indicate that they can be easily absorbed in the intestine after administration (Table 3). In the distribution section, Ellagic acid, Kievitone and Punicalin have negative values for BBB test which means they are unable to permeate through the blood-brain barrier and hence the central nervous system is protected from their action. In contrast, Remdesivir show a positive value for BBB which indicate its wider distribution and potential neurological effect. Furthermore, Ellagic acid and Punicalin were predicted to be non-substrate of P-glycoprotein which suggest both compounds could have good distribution properties. However, Kievitone and Remdesivir could act as substrate of p-glycoprotein and hence their distribution may be limited by the mediation of p-glycoprotein efflux [46].

In the metabolism part, the families of Cytp450 have been identified as an important pharmacological parameter. According to Lynch, the inhibition of these protein families may result in drug-drug interaction and bio-accumulation [36]. According to our result, the only compound found inhibiting one or more of these metabolic enzymes is Kievitone.

The excretion of all drug candidates was also assessed in our ADMET study by carrying out the biodegradation test. Interestingly, all compounds including the standard could be easily excreted as they return a negative value for biodegradation.

Finally, the pharmacodynamic profile of the selected drug candidates was also assessed. It is important to note that the failure of drugs to pass through the clinical trials is majorly because of their non-acceptable toxicological profile and hence we have put into consideration the Ames mutagenicity, Acute oral toxicity, hERG inhibition, carcinogenicity and hepatotoxicity testing in this study (Table 3). All the drug candidates and the standard drug show negative values for carcinogenicity test which may explain why the compounds may not be carcinogenic. However, it is disappointing that all the drug candidates are hepatotoxic. However, since clinically approved Remdesivir was also predicted to be hepatotoxic, it is likely that the negative consequence of this property is minimal to the liver cells. The hERG and Ames mutagenicity are crucial pharmacological parameters for evaluating if a drug-like compound could cause cardiac arrhythmia [37] and capable of mutating the DNA respectively. Ellagic acid, Kievitone and Remdesivir excel in these parameters while Punicalin with its positive values for hERG and Ames mutagenicity test may affect the rhythm of potassium channels in the heart (cardiac arrhythmia) and may induce carcinogenicity respectively. Although, Ellagic acid is the only predicted orally toxic compound in this section. However, it is evidenced that the drug candidate stands out among the selected compounds in terms of pharmacokinetic and pharmacodynamic profile. In contrast, the mutagenicity and inhibition of hERG protein by Punicalin may lead to fatal health consequence and hence we have employed a more safety measure by filtering out this phytochemical compound in our study.

Molecular dynamic is an in silico approach which is used to monitor the behavior and dynamics of binary protein-ligand complex at the atomistic level. Unlike molecular docking, MD simulation represents a more realistic delineation of the complex biological system by putting the flexibility of proteins, neutralizing ions and water environment into consideration. In this study, Ellagic acid (Ela), Kievitone (Kie) and Remdesivir (Rem) behavior with the Mpro target (6LU7) was simulated by subjecting them in Molecular dynamics environment. The reason for dropping Punicalin (the 3rd ranked inhibitor) for MD simulation is based on its poor pharmacodynamics and non-compliance with the rule of five (see Tables 3 and 2 respectively) which could be detrimental to the health of clinical trial volunteers. In order to investigate and compare the stability of 6LU7-Ela, 6LU7-Kie with the standard (6LU7-Rem), the graph spectrums of the RMSD, RMSF, H-bond and ROG of the complexes have been superimposed as depicted in Fig. 2, Fig. 3, Fig. 4, Fig. 5 .

Fig. 2

Superimposed RMSD spectrum of 6LU7 (Apoprotein), 6LU7-Ela, 6LU7-Kie and 6LU7-Rem.

Fig. 3

Superimposed RMSF spectrum of 6LU7 (Apoprotein), 6LU7-Ela, 6LU7-Kie and 6LU7-Rem.

Fig. 4

Represents the intermolecular H-bond spectrum of 6LU7-Ela, 6LU7-Kie and 6LU7-Rem.

Fig. 5

ROG Spectrum of 6LU7, 6LU7-Ela, 6LU7-Kie and 6LU7-Rem.

Root mean square deviation (RMSD)

The RMSD backbone is used to assess the structural stability of protein-ligand complex and lower value of the ligand bound protein compared to the unbound form indicates a stable system [35]. The average RMSD values for 6LU7 (Apoprotein), 6LU7-Ela, 6LU7-Kie and 6LU7-Rem are 0.2981 ± 0.0552 nm, 0.2171 ± 0.0267 nm, 0.3230 ± 0.1366and 0.2692 ± 0.0433 respectively. From this result, the RMSD value of 6LU7-Ela is lower than the apoprotein and the standard (6LU7-Rem) hence it could be regarded as the most stable complex. The spectrum of the RMSD graph also suggests minimal perturbation in the dynamical conformation of the 6LU7-Ela complex (Fig. 2). This finding may explain why Ellagic acid could not deviate away from its initial docked pose and the ruggedness to consistently bind to Mpro catalytic dyad residues. In the case of Kievitone dynamics with Mpro receptor, there is abrupt fluctuation of the RMSD spectrum between 0 and 10 ns. However, after 10 ns, stability of the complex was noticed as equilibrium was attained till the end of the 20 ns production run (Fig. 2).

Root mean square fluctuation (RMSF)

RMSF graph of the individual residues of the unbound and bound protein were also plotted to assess the flexibility of the amino acids. It is known that regions with very high RMSF values correspond mostly to the loop region which could also be involved in ligand binding and conformational alterations [33]. Interestingly, His41 (0.0787 nm) and Cys145 (0.0546 nm) catalytic dyad which are found at the loop region and believe to be critical for inhibition of Mpro were among the most stable residues of the apoprotein. Moreover, the binding of Ellagic acid, Kievitone and Remdesivir with 6LU7 return RMSF score of 0.1010 nm, 0.1141 nm and 0.0770 nm respectively for His41. The RMSF score of Cys145 for 6LU7-Ela, 6LU7-Kie and 6LU7-Rem are 0.0632 nm, 0.0583 nm and 0.0738 nm respectively (Fig. 3). Since there were no significant variations in the RMSF values of these complexes when compared with apoprotein, we conclude that the catalytic dyad residues remain stabilized. Also, the average RMSF score of 6LU7-Ela (0.1147 ± 0.0501 nm) is lower than the apoprotein, 6LU7-Kie (0.10301 ± 0.1012 nm) and 6LU7-Rem (0.1417 ± 0.0621) (see Fig. 3). Hence, it could be regarded as complex with the most stable residues.

H-bond

It is well known that intermolecular H-bond is important for maintaining the stability of protein-ligand interaction. The higher the number of H-bond formed between a complex, the stronger their interaction [41]. In this research, 6LU7-Ela averaged intermolecular H-bond formation of 1.81 for the entire 20 ns simulation run which is higher than 6LU7-Kie (1.04) and the standard (1.501) (Fig. 4). The H-bond result may explain why Ellagic acid demonstrate higher binding affinity and better stability with Mpro than both the standard (Remdesivir) and Kievitone.

Radius of gyration (ROG)

The compactness of protein in a dynamical system can be monitored by computing their ROG values. A stably folded system could be measured by equilibration of its ROG spectrum throughout the simulation period while a non-compact system may be identified when there is abrupt fluctuation of ROG spectrum [29].

The low mean ROG values of 6LU7, 6LU7-Ela, 6LU7-Kie and 6LU7-Rem (2.2267 −2544 nm) indicate minimal fluctuation of the complexes (Fig. 5). Although, the unbound 6LU7 protein is the most compact structure according to its ROG spectrum, the minimum and maximum values of 6LU7-Ela (2.21002 nm- 2.3001 nm) and 6LU7-Kie (2.20016 nm- 2.23816 nm) may explain a reasonable well folded system. Importantly, the RMSD spectrum of both complexes suggest a better compact structure when compared with the standard which shows minimum and maximum ROG values of 2.2116 nm and 2.3913 nm respectively.

The individual decomposition energy of the residues and overall binding energies of the complexes have been represented in Fig. 6 and Table 5 respectively. From the result presented in Table 5, it is crystal clear that the major driving force responsible for the spontaneity and high binding energy of the complexes is the Van der waal's interaction (VDW). In contrast, the polar solvation energy observed in the three complexes show unfavorable binding of the protein-ligand interactions. The overall result of the mm-pbsa binding energy suggests that Kievitone (−60.480±24.834 kj/mol) and Ellagic acid (−47.487±11.028 kj/mol) may be more potent inhibitor than Remdesivir(−28.911±10.9kj/mol). Finally, we plotted the graph spectrum of the individual contributing energy of the residues to assess their cruciality in the binding of the complexes. From the docked pose, His41, Met49, Cys145, Gly143, Asn142 and Glu166 were found interacting with Ellagic acid and their contributing mm-pbsa energy (in kj/mol) are 0.4035, −1.1554, −2.4237, 0.1686, 0.2228, −1.3238 respectively. Also, the contributing energy of the amino acids found interacting with Kievitone have binding energies of −0.1159 (Tyr54),2.3352 (His 41), −1.5219 (Met49), −2.1422 (Cys145), −2.796 (Met165),16.6994 (Glu166), 2.4762 (Ser144), −0.3695 (His163),0.097 (His172),−5.6508Asn142, and −0.580 (Asp187) while those binding with Remdesivir have decomposition energies of −0.6983 (His41), −0.0369 (Cys145), −4.908 (Met165), 2.9349 (Glu166), 0.1885 (Asn142), −0.7251 (His164), −0.2039 (His163), −0.202 (His172),−0.175 (Leu141), 0.055 (Gly143) and −0.079 (Ser144). From this result, it is evident that the atomic interaction of Ellagic acid and Kievitone to Cys145 is largely responsible for their significant binding energy. Interestingly, our analysis also indicates the ability of the two drug candidates to effectively inhibit SARS-COV2 main protease by binding more stably to the crucial Cys145 residue [26] than the standard.

Fig. 6

MM-PBSA Binding free energy contribution of individual residue of 6LU7 in complex with A.Ela B.Kie C.Rem.

Table 5

MMPBSA free energies of the protein-ligand complexes.

Compounds	VdW	Electrostatic	PSA	SASA	Total Binding Energy (Kj/mol)
Ellagic acid	−83.866±11.128	−34.071±7.832	78.891±9.889	−9.441±0.999	−47.487±11.028
Remdesivir	−122.849±12.780	−81.955±24.436	190.929±28.989	−15.036±1.110	−28.911±10.98
Kievitone	−126.286±40.575	−9.836±8.023	92.032±37.945	−16.390±4.363	−60.480±24.834

Quantum mechanics (DFT)

The molecular frontier orbitals are well-established parameters to provide information about the compounds' reactivity and stability. The highest occupied molecular orbital (HOMO) and the lowest occupied molecular orbital (LUMO) depict the electrons donating and receiving the compounds' ability. The energy gap, which is the difference between the LUMO and the HOMO, describes the kinetic stability and the intramolecular charge transfer. Also, compounds with significant energy gaps are referred to have low chemical reactivity and high kinetic stability. In contrast, compounds with low energy gaps are more reactive with low stability [29]. The energy gap of the lead compounds increase accordingly: Ellagic acid (0.149 a.u) < Kevitone (0.161 a.u) < Remdesivir (0.175 a.u) (Table 4 and Fig 7 ). Notwithstanding, the chemical hardness and softness were also computed from the frontier orbitals. Compounds with high chemical reactivity and low stability are considered to be soft (σ) [14], and compounds with low reactivity and high stability are referred to as hard (η) [28]. The Ellagic acid with an energy gap of 0.149 a.u has the highest chemical reactivity with less stability when compared with Kievitone and Remdesivir with an energy gap of 0.161 a.u and 0.175 a.u, respectively. While the most stable compound is Remdesivir based on the energy gap value.

Table 4

The Highest occupied molecular orbital (HOMO) energy, Lowest occupied molecular orbital (LUMO) energy.

Frontiers Orbitals	Ellagic acid (a.u)	Kievitone (a.u)	Remdesivir (a.u)
HOMO	−0.237	−0.224	−0.242
LUMO	−0.088	−0.063	−0.066
Gap (∆E)	0.149	0.161	0.175
Chemical hardness(η=(ELUMO-EHOMO)/2)	0.074	0.080	0.087
Softness (σ=1/η)	13.513	12.500	11.494

Fig. 7

The Highest occupied molecular orbital (HOMO) energy, Lowest occupied molecular orbital (LUMO) energy, and the molecular electrostatic potential (MEP).

Conclusion

In this study, we employed a computational approach with the use of molecular docking, molecular dynamics, and Quantum mechanics to investigate the inhibitory potential of dietary polyphenols against SARS-COV-2 Mpro at the molecular and atomistic levels. What we noticed from our findings is that two of the screened ligands (Ellagic acid and Kievitone) could be potent competitive inhibitors of the coronavirus main protease. It is apparent from molecular dynamics and MMPBSA analysis that the binding of Ellagic acid to His41 and Cys145 catalytic dyad of the receptor (from the docked pose) forms a more stable system when compared to the standard drug. Our findings also revealed that Kievitone is capable of disrupting the proteolytic activity of Mpro and the analysis of MMPBSA and molecular orbital theory suggest that the drug candidate is slightly stable than Ellagic acid. We, therefore, propose that Ellagic acid and Kievitone be further developed as drug candidates for the treatment of COVID-19. Nonetheless, while these findings may be appreciable and appealing under in silico research, further experimental and clinical studies are needed to buttress their preclinical and clinical relevance respectively. In addition, functional foods endowed with kievitone and ellagic acid might be functional in the management/treatment of COVID-19 disease. Although a global concern about making polyphenols functional foods is their high susceptibility towards metabolism which may impede their bioavailability index, this could be addressed by packing them into nanocarriers or through an encapsulation mechanism for efficient delivery to their target site.

Funding information

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

CRediT authorship contribution statement

Temitope Isaac Adelusi: Conceptualization, Supervision, Writing – original draft, Formal analysis. Abdul-Quddus Kehinde Oyedele: Writing – original draft, Investigation, Formal analysis, Validation. Ojo Emmanuel Monday: Writing – original draft, Methodology, Investigation. Ibrahim Damilare Boyenle: Investigation, Validation. Mukhtar Oluwaseun Idris: Software, Methodology. Abdeen Tunde Ogunlana: Methodology, Investigation, Visualization. Ashiru Mojeed Ayoola: Writing – review & editing. John Olabode Fatoki: Supervision. Oladipo Elijah Kolawole: Writing – review & editing. Kehinde Busuyi David: Validation. Akintola Adebola Olayemi: Writing – review & editing.

Declaration of Competing Interest

Authors declare no area of competing interest.