Computational studies of potential antiviral compounds from some selected Nigerian medicinal plants against SARS-CoV-2 proteins

The challenges posed by COVID-19's emergence have led to a search for its therapies. There is no cure for COVID-19 infection yet, but there is significant progress in vaccine formulation for prophylaxis and drug development (such as Paxlovid) for high-risk patients. As a contribution to the ongoing quest for solutions, this study shows potent phytocompounds identification as inhibitors of SARS-CoV-2 targets using in silico methods. We used virtual screening, molecular docking, and molecular dynamics (MD) simulations to investigate the interaction of some phytochemicals with 3CLpro, ACE2, and PLpro proteins crucial to the SARS-CoV-2 viral cycle. The predicted docking scores range from −5.5 to −9.4 kcal/mol, denoting appreciable binding of these compounds to the SARS-CoV-2 proteins and presenting a multitarget inhibition for COVID-19. Some phytocompounds interact favorably at non-active sites of the enzymes. For instance, MD simulation shows that an identified site on PLpro is stable and likely an allosteric region for inhibitor binding and modulation. These phytocompounds could be developed into effective therapy against COVID-19 and probed as potential multitarget-directed ligands and drug candidates against the SARS-CoV-2 virus. The study unveils drug repurposing, selectivity, allosteric site targeting, and multitarget-directed ligand in one piece. These concepts are three distinct approaches in the drug design and discovery pipeline.


Introduction
For over three years, the COVID-19 outbreak threatened human health across nations. Between 2002 and 2012, two highly pathogenic coronaviruses with zoonotic originssevere acute respiratory syndrome coronavirus (SARS-CoV) and Middle East respiratory syndrome coronavirus (MERS-CoV), emerged and caused devastating respiratory illnesses in humans. Consequently, a new coronavirus outbreak presented a new public health challenge in the twenty-first century. SARS-CoV-2 is the latest identified human coronavirus, the ninth of its type and the seventh discovered within the past two decades [1]. Although the origin of SARS-CoV-2 has not been fully established, it was first reported (in late December 2019) after several patients showed pneumonia of unknown origin from multiple healthcare facilities in Wuhan, Hubei Province of China [2].
The primary rapid spread of the virus is through transmission via droplets from the lung parenchyma of infected persons during coughing or sneezing [3]. The spike protein of the COVID-19 virus interacts with the host cell enzyme furinavailable in several tissues, thus giving the virus capacity to propagate faster and affect several organs [4]. The COVID-19 virus infection exhibits several symptoms through the respiratory system ensemble of ear, nose, and throat (ENT). COVID-19 infection further spreads to attack cardiovascular, gastrointestinal, musculoskeletal, ocular, cutaneous, neurological, and hematological systems. So, symptoms include pneumonia, acute respiratory distress syndrome (ARDS), and multi-organ dysfunction [5]. Therapeutic options for COVID-19 treatment include antivirals, antiparasitic, anti-inflammatories, interferon, and convalescent plasma. Monoclonal antibodies, hyper immunoglobulin RNAi, mesenchymal stem cell therapy, and vaccines are also therapeutic options [6]. Aside from vaccines, which cannot suffice in treating already infected individuals, the U.S. Food and Drug Administration (FDA) has approved only one pharmaceutical product for COVID-19 infection treatment. As of the inception and completion of this study, the were no approved anti-COVID molecule. So, no information on the approved drug is available or discussed in this study.
Studies have shown significant reduction/disappearance of viral load in COVID-19 patients treated with hydroxychloroquine and azithromycin [7]. Similarly, lopinavir-ritonavir, favipiravir, and remdesivir are suggestive COVID-19 drugs. However, on 04 July 2020, the World Health Organization (WHO) ruled on using hydroxychloroquine and lopinavir-ritonavir as anti-COVID-19 drugs following solidarity trials of the suggested therapies [8]. The traditional method of developing new drug entities requires a thorough scientific evaluation of potency, efficacy, bioavailability, adverse effect on non-target sites, safety, and preclinical and clinical tests [5,9]. This process requires substantial funds and usually takes months/years before a single drug is certified for clinical use [10], thus making drug repurposing and screening for therapeutic options from medicinal plants a viable option. Medicinal plants have served as the basis for traditional medicines in treating various diseases in diverse cultures of the world. Its use is as old as human civilization [11]. Nigeria has a rich biodiversity of plants with medicinal value. Several plants from Nigeria have been identified to possess antiviral properties and are commonly used in traditional medicines [12]. Therefore, it would be logical to evaluate these plants from which enormous research findings and traditionally used antiviral properties exist.
The discovery and development of potent drug leads are plausible by applying computational approaches like machine learning and docking [13]. Scientists have used molecular docking to study the interaction of proteins and drug molecules. According to Gao et al. (2020), machine learning and molecular docking techniques create added value to data reported by other scientists [14]. Due to the reuse of free access data and available software, we can use computational methods to validate, test, and provide a hypothesis, thereby considerably reducing research costs [15]. With the available data, this research features the interaction of active compounds in known medicinal plants against three established targets of SARS-CoV-2. These proteins are 3-chymotrypsin-like protease or main protease, human angiotensin-converting enzyme 2, and papain-like protease. The proteins hereafter referred to as 3CLpro, ACE2, and PLpro perform various functions in the SARS-CoV-2 virus cycle [16]. For instance, 3CLpro cleaves the SARS-CoV-2 virus polyprotein at eleven distinct regions [17], ACE2 speeds the catalysis of angiotensin I to angiotensin-(1-9) [18], and PLpro processes the viral polyproteins enabling viral spread [19].

Plant and compound selection for ADME screening
The selection of the four plants (Afromomum meleguata, Vernonia amygdalina, Cajanus cajan, and Arctostaphylos uva-ursi) premised on their ethnopharmacological history from the Nigerian medicinal database (http://www.medicinalplantsinnigeria.com). Also, we are interested in those with antiviral properties guided by previous investigations. The search reveals fifty-four isolated phytocompounds from these plants. We screened these compounds for pharmacokinetics properties using an integrative in silico model. These properties include absorption, distribution, metabolism, and excretion (ADME) to screen for natural compounds that may be bioactive via oral administration.

Ligand and receptor preparation for molecular docking
Of the 54 compounds, 30 showed high gastrointestinal (GI) absorption, 15 are potent antiviral candidates, and more than 5 are potential anti-COVID-19 compounds. We docked the 15 phytocompounds ( Fig. 1) with the most promising ADME and drug-like properties into the active sites of 3CLpro, ACE2, and PLpro of the SARS-CoV-2 virus. The process involves energetically minimizing the retrieved ligand structure using the AutoDock Vina tool [37] to obtain the PDBQT format. The three-dimensional structures of SARS-CoV-2 Main protease, ACE2, and PLpro are available in the Protein Data Bank (PDB) [38] with codes 6LU7 [39], 7DMU [40], and 7LBS [41]. We remove all water molecules and heteroatoms from the crystal structures using Discovery Studio [42] and save them in PDB format for further analysis.
AutoDock [43] facilitates protein optimization, polar hydrogen group addition, and inducing charges to the ligand and protein structure atoms through AutoDockTools. The grid box size for the selected phytocompounds interaction within the various protein entails mapping out the active sites using the experimental information on the correct active sites for ligand binding. We used AutoDock to couple these phytochemicals to the active sites of 3CLpro, ACE2, and PLpro, and selected the accurate docking approach. The population size (150), maximum energy evaluation (2,500,000), number of generations (27,000), and other terms follow the default Lamarckian genetic parameterization of the AutoDock [43] software. The output shows the ten most favorable binding poses and the reported outcome is the most suitable for each simulation.
To explore if these compounds attach favorably at the active sites or perhaps allosteric sites on these targets, we used SwissDock [44,45] to predict the poses of these compounds on the proteins in an unconstrained but accurate docking procedure. Using more than one docking software enabled us to ascertain that the outcome is not a bias or artifact of one docking program. We reported the best docking score in each case and viewed the interaction modes in the Discovery [42] software. The SwissDock algorithm integrates several other programs, facilitating accurate molecular docking within minutes [44,45]. The software can yield robust binding modes up to 15,000. The reported energy per each protein-ligand docking simulation is the most favorable outcome from the clusters dumped into the result file after energy ranking with solvent consideration using the FACTS [46] implicit solvation scheme and clustering.

Molecular dynamics simulations
We performed MD simulations in the AMBER 18 program package with an integrated graphics processing unit (GPU) [47] using the particle mesh Ewald method and the ff19SB force field [48] for the protein.
The periodic boundary condition has long-range electrostatic interactions [49] at a 12 Å cut-off. We generated a partial atomic charge for the ligand using the AM1BCC [50] scheme of the ANTECHAMBER. The AM1BCC is suitable for high-quality charge assignment on small molecules [50]. Further steps include topology and coordinate preparation, enzyme-inhibitor coupling, neutralization, and solvation. We explicitly solvated the complex in the TIP3P orthorhombic water box at 10 Å to any edge using the Leap module of the AMBER 18 package.
We used 10000 steepest descent steps with 10 kcal/mol/Å 2 harmonic restraint on all heavy atoms to relax the system and remove potential atom clashes. The complete minimization entails another steepest and conjugate gradient descent of 5000 steps without restraining. The enzyme-ligand was heated gradually from 0 K to 300 K for 300 ps using Langevin dynamics with 1ps collision frequency, and a 5 kcal/mol/Å 2 applied harmonic restrictions at a constant volume. We subsequently equilibrated for 500 ps at 300 K under a fixed pressure and temperature (NPT) ensemble. The final step was MD simulation for 120 ns at 2.0 fs time step with no restraints at 300 K and 1 atm in the NPT ensemble using the Langevin temperature scaling [51] and Berendsen barostat [52] to maintain the temperature and pressure, respectively. The production simulation also entails applying the SHAKE algorithm [53] to constrain all bonds involving hydrogen atoms.
We measured system stability through root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), the radius of gyration (RoG), and solvent-accessible surface area (SASA) calculations. RMSD, RMSF, and RoG calculations are from equation (Eq.) 1-3 integration in the CPPTRAJ [54] module implemented on the AMBER package. The RMSD of the protein backbone alpha carbon (Cα) represents the standard deviation (Eq. (1)) of the interatomic distance between Cα backbone atoms of two amino acids [55]. The v i in Eq. (1) signifies the coordinates of the Cα atom in v at the time i, and w i denotes the coordinates of the Cα atom in w at the time i.
RMSF is like RMSD, but on a per-residue basis, it presents a unique way to predict the conformational changes of the protein scaffold. RMSD facilitates protein motion prediction over the simulation time, while RMSF measures flexibility per atom or residue for all trajectories. In Eq.
(2), x i(j) represents the ith atom position in the jth residue, and (x i ) denotes the averaged location of the ith atom in all residues.
The RoG is the moment inertia of atoms from their center of mass, which quantifies the molecular rigidity [56]. The RoG (Eq. (3)) is the square root of the inertia moment of atoms where n is the atom number, r i represents the atomic position, and r m indicates the mean position of all the elements.
The Molecular Mechanics/Generalized Born Surface Area (MM/ GBSA) method of the MMPBSA.py program implemented in AMBER 18 was used to estimate the relative binding free energies for the dynamical enzyme-inhibitor system. Eq. (4) illustrates how the MM/PBSA.py program computes the binding free energy (ΔG bind ). ΔE MM and ΔG sol are the molecular mechanical interaction and the solvation-free energies, respectively, TΔS denotes the entropy contribution. The sum of the van der Waals (ΔE vdw ) interaction and the electrostatic (ΔE ele ) interaction energies comprise the molecular mechanics' interaction energy (ΔE MM ) between the receptor and ligand (Eq. (5)). The relative binding energies were calculated using 1200 frames from the 120 ns, with the conformational entropy via the AMBER NMODE module [57,58]. The input file for the MM/GBSA calculation loops over the entire 120,000 trajectories by skipping every 100 frames to output the binding energy components.
The nonpolar term of the solvation energy was determined through the SASA (Eq. (7)), where β and γ are set at − 0.5692 kcal/mol and 0.0378 kcal/mol/Å 2 , respectively. The implemented water probe radius for SASA determination is 1.4 Å. Further analysis includes linear interaction energy (LIE) of each ligand with the entire protein scaffold over the 120 ns simulation.

ADME analysis
In this work, we select four plants based on their phytocompound bioavailability for virtual screening. The plants Afromomum meleguata, Vernonia amygdalina, Cajanus cajan, and Arctostaphylos uva-ursi have various bioactive components, including antioxidant, antiviral, and antimalarial properties when consumed orally. The screening process of the plant thus requires that the available antiviral moiety is absorbable via oral prescription. Therefore, the herbs should contain biologically proven antiviral phytocompound with sufficient drug-likeness and oral bioavailability evaluations.
There are 54 bioactive compounds present in the plant after crosschecking with the Dr. Duke plant database (https://phytochem.nal. usda.gov/phytochem/search). This analysis enables us to identify the number of phytocompound present in the plant. We screen these phytocompounds on SwissADME to estimate absorbability, distribution, metabolism, and discharge (Table 1). These properties are necessary recipes in drug discovery or pharmacokinetics [59]. Among the 54 phytocompounds, 30 showed high GI absorption but had no implications as antiviral agents, and 15 passed the screening test showing high ADME. A literature survey shows that the identified compounds have antiviral components capable of inhibiting SARS-CoV-2 or other contagious respiratory ailments.

Molecular docking of phytochemicals in SARS-CoV-2 virus proteins
After obtaining the phytocompounds with antiviral potentials through screening, we subjected them to a molecular docking process to reduce the risk of blind ligand prediction devoid of target interaction. The predicted best binding scores using AutoDock/SwissDock are available in Table 2, the 15 phytocompounds showing favorable interactions with the studied SARS-CoV-2 virus targets. The binding energy is within − 5.5 and − 9.4 kcal/mol through the docking programs, indicating these phytocompounds as potent inhibitors of 3CLpro, ACE2, and PLpro targets of the SARS-CoV-2 virus. The unconstrained accurate docking model with SwissDock shows some of these compounds bind favorably at other sites than the active sites of the proteins (Fig. 2). To verify the accuracy of the docking procedure, we re-docked the cocrystallized inhibitors of all the targets using SwissDock. An overlay of the most favorable docked pose of the inhibitors with the retrieved x-ray model shows that the blind docking detects the active sites, and the other sites are potent allosteric sites (Fig. S1 in the Supporting Information, SI). Zoomed-in interaction profiles of the phytocompounds with several regions of the 3CLpro, ACE2, and PLpro targets of the SARS-CoV-2 virus are available in Fig. 3, 4, and 5. Each snapshot of the nonbonded residue interaction with the ligand is at a 4 Å neighbor distance.
Although these compounds bind at the active sites in all cases, some attached favorably at distal sites away from the active sites (Fig. 2). For instance, isorhamnetin binds better at allosteric sites than the active regions in 3CLpro and PLpro (Figs. 3 and 4). The estimated binding scores are − 7.8 and − 7.4 kcal/mol, respectively ( Table 2). The corresponding interaction energy values for isorhamnetin coupling at the active sites of 3CLpro and PLpro are − 7.3 and − 7.0 kcal/mol ( Table 2). Isorhamnetin interacts best at the ACE2 active site with binding affinities of − 8.2 and − 7.8 kcal/mol. Various studies have provided the scientific basis for the health benefit of isorhamnetin. Isorhamnetin is a flavonoidavailable in eukaryotic cells and plants like parsley, green bell peppers, and dills. It is also available in a lower concentration in Romaine lettuces, Chinese cabbages, and pears [60]. Isorhamnetin occurs in various forms, with each derivative having unique bioactivity. For example, isorhamnetin-3-O-glucosyl-rhamnoside exhibits significant anti-inflammatory effects [61]. In a study on several substituted isorhamnetin, researchers [62] noted appreciable binding stability of these moieties to the SARS-CoV-2 3CLpro. The derivatives bind better than the standard drug darunavir. Our outcome indicates that isorhamnetin, with moderate solubility in water and overall appreciable drug-like properties (Table 1), is a potent compound that can prevent different steps in the SARS-CoV-2 virus life cycle, thus having potential antiviral effects against SARS-CoV-2.
All 15 phytocompounds bind appreciably at ACE2 active site (Fig. 2) with a very high binding affinity (Table 2). It is rational to suggest that these compounds show more selectivity for ACE2 than 3CLpro and PLpro. On average, both docking protocols show favorable interaction between ACE2 and the phytocompounds. 2-hydroxyl genistein shows a better affinity for ACE2 than genistein, with an almost 2 kcal/mol difference in the constrained docking (Table 2). Hydroxyl genistein perturbs bioactivity significantly with more favored Lipinski's rule of five, ADME, and drug-like properties than its analog (Table 1). Meanwhile, there is an outcome of a computational study showing genistein as a potent inhibitor of SARS-CoV-2 main protease with a docking score of − 7.6 kcal/mol [26]. This value is closer to − 7.4 kcal/mol obtained with SwissDock for this molecule interaction at the active site of ACE2 (Table 2).

Table 2
Docking scores (binding energy in kcal/mol) of the phytocompounds against SARS-CoV-2 proteins using AutoDock/SwissDock (as the values appear) for the best binding mode. The words activeS and allos.S imply active site and allosteric site, respectively, for the most favorable binding energy per docking. These items are not indicated for compounds with favorable active site binding in both docking programs. phytocompounds (Fig. 2) at the proposed ACE2 active site [65] aligns with the feasibility of this region as the appropriate catalytic region of ACE2. Despite similarities, the studied phytocompounds exhibit some differing characteristics. For instance, apigenin and emodin are structurally related (Fig. 1) with the same molecular weight (Table 1) but with different binding potentials. Apigenin binds at ACE2 active site with − 8.1 kcal/mol affinity, whereas emodin shows − 9.2 kcal/mol for the same interaction (Table 2). Emodin is an anthraquinone moiety available in several plants' roots, leaves, and bark, including aloe vera, cascara, rhubarb, and senna. In traditional medicine, emodin has many beneficial effects like a diuretic, vasorelaxant, anti-bacterial, antiviral, anti-ulcerogenic, anti-inflammatory, and anti-cancer. These properties enable various disease management, including cardiovascular diseases and osteoporosis [66]. Herein, the derived emodin significantly binds the selected SARS-CoV-2 targets, so it could potentially block spike protein interaction with ACE2. The approximate binding energy of emodin coupling with the ACE2 active site is − 9.2 kcal/mol (Table 1).
Apigenin is a natural flavonoid (4 ′ ,5,7,-trihydroxy flavone) found in fruits (oranges), vegetables (onions, celery), spices (basil, oregano), and herbs (tea, beer, and red wine). It is in class II per the biopharmaceutics classification system. Besides its identification as a potent 3CLpro inhibitor at the molecular level [22,67], our screening identifies apigenin with favorable binding to all the studied targets (Table 2). Apigenin, 3-O-methyl quercetin, boeravinone H, isorhamnetin, luteolin, and quercetin show approximately the same affinity (− 7.0 kcal/mol) to PLpro active site with AutoDock. The outcome of the unrestraint docking is fascinating. The overall interaction energies of these phytocompounds to the PLpro active site using AutoDock and SwissDock are within − 5.6 to − 8.0 kcal/mol ( Table 2). Only five compounds (boeravinone H, cajaisoflavone, capsaicin, luteolin, and salidroside) converge consistently at the active site. Eight phytochemicals maintain stability at the N-terminal ubiquitin-like (Ubl) domain, while two others (2-hydroxyl genistein and formononetin) bind around the catalytic domain (Fig. 2). Earlier experimental studies on the MERS PLpro elucidation showed the possibility of the Ubl domain binding molecules appreciably [68].
Although salidroside failed Ghose drug predicts theory, it has a laudable outcome. Salidroside is the most soluble and permeable (Table 1); it binds at PLpro active site with − 5.8/− 7.8 kcal/mol for the AutoDock/SwissDock screening ( Table 2). Fig. 5 shows the interaction modes for the five possible binding regions identified using blind molecular docking. Emodin is among the seven compounds that favorably bind at an established allosteric site of PLpro through the Ubl domain [68], while genistein binds preferably at another region of this domain (Fig. 5). Emodin and genistein show binding affinities of − 7.1 and − 7.5 kcal/mol, respectively ( Table 2). We attribute this difference to the many possible HB interactions in genistein compared to emodin, with only one hydrogen bonding through Arg65 of PLpro (Fig. 5). Researchers [69] have proposed emodin as MTDL of anti-COVID-19 capable of interacting with several pathways, including ACE2 and 3CLpro.
The boeravinone class, including A-J, are potent inhibitors of cancer targets [76]. There are suggestions that boeravinone H is a substantial component of Boerhavia diffusa with the potential to mop up the viral particles, thereby preventing the initial phase of the hepatitis C virus attachment [77][78][79]. Review studies have featured Boerhavia diffusa with boeravinone active components as a potential candidate for COVID-19 management [80] and RNA virus [81] inhibition. Mitra et al. (2020) identified boeravinone F in an in silico study as a potent 3CLpro and PLpro dual inhibitor [82] in the SARS-CoV-2 virus. Therefore, boeravinone H identification in our screening and docking approach looks promising as a likely anti-COVID-19 molecule. The identified boeravinone H (Fig. 1) has noteworthy ADME, drug-like, and physicochemical properties ( Table 1) that favor its significant interaction with 3CLpro, ACE2, and PLpro targets of SARS-CoV-2 virus within − 6.9 and − 9.3 kcal/mol binding energy (Table 2).

Probing the allosteric sites of SARS-CoV-2 virus proteins
Since there are several predictions on the binding sites of 3CLpro, ACE2, and PLpro, we explore the stability of the identified non-catalytic binding sites of these proteins that are unavailable in the literature. To further evaluate if the phytochemicals bind preferably at these regions over the active site, we simulate the most favored compound-allosteric site interaction for the catalytic site using the dynamics of the uninhibited protein as a control. So, we have a total of 1.440 μs simulations for 3CLpro (3 systems), ACE2 (4 complexes), and PLpro (5 systems). Table 3 shows the binding energy profile of inhibitor-bound protein structures, whereby we also consider only the most favored potent allosteric compound at the active site.
The crucial energy component behind the favorable non-bonded interaction of these compounds to the SARS-CoV-2 virus proteins is the vdW energy (Table 3), including the various aromatic and hydrophobic interactions (Figs. 3-5). Cajaisoflavone-ACE2 allosteric interaction shows the highest (− 45.8 kcal/mol) contribution, while genistein contact at PLpro active site shows the lowest vdW contribution of − 25.25 kcal/mol. HB interactions are signatures of the electrostatic energy driving potent molecules to the targets [49]. Cajaisoflavone interaction at ACE2 active site shows the most favored ΔE ele outcome of − 15.5 kcal/mol, while genistein-PLpro active site ΔE ele is the least favorable with − 5.7 kcal/mol. The average binding energy of the systems in the gas phase ranges from − 31 to − 60.5 kcal/mol with the MM/GBSA scheme (Table 3), while the outcomes with LIE code are within − 25.7 and − 60.4 kcal/mol (Fig. 6). Solvent contribution to the binding averaged within 11.3 and 18.8 kcal/mol, which perturbs the overall binding energy (ΔH). Cajaisoflavone interaction at the active site of ACE2 shows the most favored ΔH and TΔS ( Table 3).
The probability distribution of the total interaction energy for the protein-ligand binding shows that the selected phytocompounds maintained bound states throughout the simulation, ranging from 0 to − 80 kcal/mol (Fig. 6). Formononetin-PLpro allosteric site system has a unique bimodal distribution, isorhamnetin-3CLpro (active site), cajaisoflavone-ACE2 (active region), and boeravinone H-ACE2 (nonactive site) are perfectly unimodal, while others have pseudo unimodal distribution. It is intriguing to notice that the interactions of these compounds at the various protein active sites have relatively significant values with the tallest peaks and narrowed densities depicted with black, blue, and yellow legends (Fig. 6) for isorhamnetin-3CLpro, cajaisoflavone-ACE2, and genistein-PLpro systems, respectively. The energy values (in kcal/mol) − 40.5, − 60.4, and − 32.7 for the respective protein active binding show the selectivity of these compounds for the intrinsic active region. Only PLpro allosteric interaction with genistein at a non-catalytic site shows more favored interaction energy (− 40.4 kcal/mol) than the active site counterpart (− 32.7 kcal/mol). All the ACE2 systems have the highest average interaction energy values. Overall, the average coupling energy trend for the studied ligands from LIE calculations (Fig. 6) matches the MM/GBSA (ΔG gas ) analysis (Table 3). Cajaisoflavone at ACE2 active site and 2-hydroxyl genistein-PLpro non-active site have the highest and lowest values, respectively, in both energy analyses.
All the complexes and apoenzymes showed unique behavior over the simulation time (Fig. 7). The alpha carbon RMSD plots relative to the starting PDB structure show that all the proteins' primary structures attain various conformational states during the simulation. Fig. 7A depicts the 3CLpro Cα atoms of the protein apo structure, active site, and non-catalytic bound states with isorhamnetin, showing that the protein backbone has at least two conformational states (from the histogram distribution) during the simulation. The average RMSD of the 120 ns simulated SARS-CoV-2 virus main protease is approximately 2 Å relative to the initial PDB structure for the three possibilities (Fig. 7A). Besides the allosteric site outcome with a separation midway through the simulation, the RMSD values overlap, denoting similarity in the conformational evolution of the 3CLpro backbone per time (Fig. 7A). Analysis of the relative ligand RMSD of all atoms is available in Fig. S2 of the SI. The output shows that the phytochemicals are stable within their bounded region throughout the simulations compared to the input conformation.
There is no wide separation in the rigidity, measured with RoG, of the main protease at unliganded or bound states (Fig. 7B). In contrast, the apo 3CLpro has the highest solvent accessible surface area of 14,616.6 Å 2 , indicating that the solvent has more surface to occupy than the bound states. The calculated SASA values (Fig. 7C) for isorhamnetin active and non-active site 3CLpro systems are 14,034 and 14,193.5 Å 2 , respectively. 3CLpro complexes show distinct flexibility on the perresidue scale (Fig. 7D) with average values approximated to 1.4, 8.7, and 26.8 Å for apo, isorhamnetin active, and non-active, respectively. Isorhamnetin is the only molecule that binds favorably at a non-active  ΔE vdW = van der Waals energy; ΔE ele = Electrostatic energy; ΔG gas = Gas phase energy; ΔG solv = Solvation energy; ΔH ¼ total binding energy or enthalpy; TΔS = entropy contribution; ΔG bind ¼ binding free energy; ActiveS = active site; and allos.S = allosteric site.
site of 3CLpro (Fig. 3), but its binding perturbs 3CLpro flexibility in 25folds (Fig. 7D). However, isorhamnetin binding to the 3CLpro active site is more significant than at the non-active site with about − 3.5 kcal/mol binding free energy difference (ΔG bind , Table 3). Based on the outcomes, isorhamnetin in 3CLpro is more suited as an active site molecule than an allosteric binder. However, this observation does not outweigh its probability as a potential allosteric 3CLpro inhibitor. The alpha carbon dynamics in ACE2 systems (Fig. 7E) have almost the same outcome as 3CLpro RMSD (Fig. 7A) with over two conformational states, except cajaisoflavone-active site interaction showing a unimodal distribution. The protein backbone shows similar incremental dynamical projections per time (Fig. 7E). Cajaisoflavone at ACE2 active site maintained a stable Cα atom evolution over the simulation relative to the initial structure and has the lowest deviation averaged to 1.75 Å. The highest average RMSD (2.55 Å) in apo ACE2 denotes a large-scale motion in the protein backbone over the simulation period. Ligand coupling to its active site lowers this projection significantly, while nonactive site inhibition slightly decreases the Cα atom flexibility. Using RMSD as a stability metric, it is logical to conclude that cajaisoflavone at ACE2 active site maintains relative stability throughout the simulation.
The RoG and SASA projections in ACE2 ( Fig. 7F and G) are distinct with increasing values as the simulations proceed, unlike the 3CLpro and PLpro systems with almost even distribution. However, the average RoG for ACE2 systems is close (Fig. 7F), while the solvent-accessible area within the boeravinone H non-active site complex is higher than others (Fig. 7G). Interestingly, the flexibility of the protein on the per-residue scale (Fig. 7H) differs from its backbone analysis (Fig. 7E). The apo ACE2 shows a relatively lower RMSF (10.92 Å), and cajaisoflavone at its active site yielded the lowest average value of 10.81 Å, denoting that ligand at ACE2 active site restricted the protein's mobility. In contrast, the ACE2 residues are highly flexible for the allosteric systems, with average RMSF scores of 19.4 and 20.5 Å (Fig. 7H) for cajaisoflavone and boeravinone H, respectively. The calculated binding free energy of the phytocompounds to ACE2 active and non-active sites are favorable with values − 11.5, − 5.1, and − 1.5 (Table 3). Besides the dynamical metrics RMSD, RoG, and RMSF showing appreciable outcomes for ACE2cajaisoflavone (active site), several energy metrics, including docking score (Table 2), ΔG gas , ΔH, TΔS, and ΔG bind (Table 3) favor the interaction.
The trajectories of the PLpro systems have almost uniform projections in RMSD, RoG, and SASA ( Fig. 7I-K) analyses. However, a closer look shows alpha carbon dynamics attains two distinct conformational states when formononetin binds at PLpro (non-active site) compared to the other four systems with one bound state (Fig. 7I). The Formononetin system also yields the highest average RMSD value of 2.07 Å, while 2-hydroxyl genistein shows the lowest (1.69 Å) outcome (Fig. 7I). The evolution of residue rigidity per time, analyzed with the radius of gyration, depicts almost the same projection for the PLpro systems with the lowest value of 23.65 Å from genistein non-catalytic site interaction (Fig. 7J). Similarly, the available surface within the PLpro for solvent is the lowest (14,632.4 Å 2 ) in the genistein-allosteric site bound (Fig. 7K). This observation is almost the same for the RMSF analysis, where genistein-PLpro (non-catalytic site) interaction shows the lowest average (~16 Å) per-residue flexibility among other bound states. PLpro apo structure maintained the least mobility (Fig. 7K) like the 3CLpro counterpart (Fig. 7D).
The induced residue dynamics by these compounds vary per the region of interaction on the protein. For instance, formononetin interaction at a non-catalytic area on the PLpro perturbs flexibility around residues 1 through 75 (Fig. 7L) because it binds far away from this domain (Fig. 5). The MM/GBSA and linear interaction energy calculations favor genistein (allosteric site) over other experiments with Fig. 6. Total interaction energy of phytocompounds bound states at the active and allosteric sites of SARS-CoV-2 virus proteins 3CLpro, ACE2, and PLpro for each 120 ns production run.
significant ΔG gas , ΔH, TΔS, ΔG bind (Table 3), and interaction energy (Fig. 6). The total binding free energy (PLpro-genistein, allosteric site) is approximately − 3 kcal/mol, but its interaction at the catalytic site shows ~8 kcal/mol. These results buttressed the lower active site-constrained docking score of − 5.6 kcal/mol than the non-active site with − 7.5 kcal/mol ( Table 2). Based on the studied metrics, genistein is relatively stable with several pointers to be a potential allosteric inhibitor of PLpro. There are emerging possibilities to establish new methods to recover these phytocompounds from alternative sources like macro-and microalgae [83,84], plus agro-food waste routes [85]. These green and innovative approaches, including enzymatic treatment and microwave-assisted extraction, are environmentally friendly and sustainable. Also, ultrasound-assisted extraction, supercritical CO 2 , and subcritical water technologies are intriguing green methodologies [86,87]. All the identified phytocompounds fit Lipinski's and other drug-likeness rules and fulfill the necessary properties as potent bioactive molecules. Therefore, we recommend their development into effective therapy against COVID-19 and further exploration as a potential drug candidate against the SARS-CoV-2 virus.

Conclusion
For centuries, medicinal plants and foods remained critical sources of efficient and safe consumables for humans and animals. Among them are flavonoids, widely distributed in vegetables, fruits, and leaves. They are known for their anti-inflammatory, antibacterial, and anti-tumor effects [88,89]. They are cheap, safe, and effective inhibitors and could help as effective compounds against the COVID-19 pandemic. Using computational and theoretical methods, we have identified some potent antiviral inhibitors that can serve as a starting point to stop the replication of SARS-CoV-2. These naturally occurring compounds derived from known Nigerian medicinal plants show appreciable affinities for 3CLpro, ACE2, and PLpro of SARS-CoV-2. The virtual screening and molecular docking procedures enabled us to elucidate some crucial drug discovery paradigms on one footing. The identified compounds have antiviral implications, thereby satisfying drug repurposing. The molecules are potent multitarget inhibitors of three key targets of the SARS-CoV-2 virus. These phytocompounds are selective for the active site, while some show potential as allosteric inhibitors, including genistein. Hence, this study examines drug repurposing, selectivity, MTDL design, and allosteric binding. Among the 54 phytocompounds screened from the four medicinal plants based on their ADME and antiviral properties, 15 showed good antiviral properties. Improved analysis with the MD study showed that these phytochemicals bind against 3CLpro, ACE2, and PLpro of the SARS-CoV-2 virus. The identified compounds include the popular naringenin and quercetin, with appreciable binding energy to the selected targets. Some compounds like cajaisoflavone and isorhamnetin are selective for 3CLpro and ACE2 active sites, respectively, while genistein favorably binds at the non-catalytic domain of PLpro. Our study shows that virtual screening, molecular docking, and molecular dynamics simulations are critical to drug prediction and bioactive molecule identification as potential drug candidate.

Funding
MML is grateful for the financial support of the National Research Foundation of South Africa (Grant Number: 120707).

Ethical approval
Not applicable.