In Silico Evaluation of Different Flavonoids from Medicinal Plants for Their Potency against SARS-CoV-2

The ongoing pandemic situation of COVID-19 caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) poses a global threat to both the world economy and public health. Therefore, there is an urgent need to discover effective vaccines or drugs to fight against this virus. The flavonoids and their medicinal plant sources have already exhibited various biological effects, including antiviral, anti-inflammatory, antioxidant, etc. This study was designed to evaluate different flavonoids from medicinal plants as potential inhibitors against the spike protein (Sp) and main protease (Mpro) of SARS-CoV-2 using various computational approaches such as molecular docking, molecular dynamics. The binding affinity and inhibitory effects of all studied flavonoids were discussed and compared with some antiviral drugs that are currently being used in COVID-19 treatment namely favipiravir, lopinavir, and hydroxychloroquine, respectively. Among all studies flavonoids and proposed antiviral drugs, luteolin and mundulinol exhibited the highest binding affinity toward Mpro and Sp. Drug-likeness and ADMET studies revealed that the chosen flavonoids are safe and non-toxic. One hundred ns-MD simulations were implemented for luteolin-Mpro, mundulinol-Mpro, luteolin-Sp, and mundulinol-Sp complexes and the results revealed strong stability of these flavonoid-protein complexes. Furthermore, MM/PBSA confirms the stability of luteolin and mundulinol interactions within the active sites of this protein. In conclusion, our findings reveal that the promising activity of luteolin and mundulinol as inhibitors against COVID-19 via inhibiting the spike protein and major protease of SARS CoV-2, and we urge further research to achieve the clinical significance of our proposed molecular-based efficacy.


Introduction
COVID-19, manifested by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has found its place in the concept of a new emerging outbreak that was first recorded in patients suspected with symptoms of the possible common cold found in Wuhan city of China in late December 2019, followed by a rapid spread worldwide until now [1]. SARS-CoV-2 from the coronaviridae family possesses a single-stranded RNA molecule that infects both humans and mammals. The clinical manifestation covers a wide range of symptoms including lymphopenia, fever, dry cough, diarrhea, rhinorrhea, vomiting, nausea, and fatigue [2,3]. The SARS-CoV-2 is comprised of four main viral proteins such as spike glycoprotein (S), envelope protein (E), membrane protein (M), and nucleocapsid protein (N) [4]. The human angiotensin-converting enzyme 2 (HACE2) is a receptor on the host cells vital for the incorporation of SARS-CoV-2 spike protein into the lung tissue [5,6]. Many studies have revealed the crystal structures of the main protease of SARS CoV-2, including or excluding inhibitory molecules. Thus, it has been the most characterized target for novel drug discovery [7][8][9]. To design novel drugs against SARS CoV-2, we have been provided with several efforts to obtain an understanding of the protein-ligand interaction on a molecular basis [10]. Currently, a potential therapeutic drug is yet to be discovered; however, it is reported that chloroquine, hydroxylchloroquine, favipiravir, lopinavir, simeprevir, darunavir, and remdesivir can act as a helpful medication against this virus [11][12][13]. Additionally, the combination of lopinavir and ritonavir with nebulized alfainterferon is described as a sympathetic medication against SARS-CoV-2 [14][15][16][17]. The role of traditional medicine in COVID-19 treatment is studied and reported [18][19][20]. Medicinal plants containing bioactive compounds have always been of great interest to scientists as they play a significant role in curing human diseases. Therefore, they can reinforce the existing drugs to elicit improved effects against different diseases without negligible adverse effects [21]. Over 80% of the African population relies on traditional medicines against different pathogens such as primary health problems according to WHO [22]. Flavonoids originate from medicinal plants, such as phenolic compounds, and have several biological applications including antiviral, anti-inflammatory, antibacterial, antifungal, anti-aging, anti-Alzheimer's disease, antitumor, antioxidant, etc. [23][24][25][26][27][28][29][30][31]. In addition, flavonoids exhibit antiviral activity against various viruses such as Sindbis virus, HCMV, HSV-1, HSV-2, DENV, parainfluenzavirus 3, chikungunya virus [32][33][34][35][36][37][38]. Among the naturally occurring flavonoids, anthranol, apigenin, cannflavin, derrisin, jaceidin, lupiwighteone, luteolin, methylglovanon, mundulinol, naringenin, rhamnetin, and tamarixetin, and due to their significant biological roles they are reported the most important flavonoids [39][40][41][42]. These flavonoids are abundant in various medicinal plants such as Andrographis paniculata, Moringa peregrina, Santolina insularis, Melissa officinalis, Solanum nigrum, and Tamarix nilotica and these plants are already used as traditional medicine in some countries [43][44][45][46][47][48]. Therefore, the present study aimed to depict the inhibitory potential of these flavonoids with variable pharmacological effects, targeting the spike protein (Sp) and main protease (M pro ) of SARS-CoV-2 by in silico approaches. The inhibitory effects of all studied flavonoids were discussed and compared with some currently used antiviral drugs in COVID-19 treatment.

Molecular Docking
The binding interaction obtained due to the interactions among all the studied compounds inside 6W63 and 6VYB has been displayed in Figures 1-4, respectively. Table 1 shows the result obtained from the molecular docking analysis for three proposed antiviral drugs (Favipiravir, Lopinavir, and Hydroxychloroquine) and the studied compounds against 6W63. The molecular docking analysis revealed that only luteolin and mundulinol have the highest binding energies (−13.97 and −13.12 kcal/mol, respectively) compared to favipiravir, lopinavir, and hydroxychloroquine in the state of 6W63 (−13.02, −12.98, and −12.45 kcal/mol, respectively). Different types of interactions such as hydrogen bond, Van der Waals, and pi-alkyl interactions were established between the amino acid residuals of 6W63 and the studied compounds with the highest binding energy.  PRO 52, and TYR 54 amino acid residues through Van der Walls interactions, whereas in the state of 6VYB, the molecular docking analysis showed that mundulinol, derrisin, and luteolin compounds displayed the highest binding energies (−11.08, −11.04, and −10.92 kcal/mol, respectively) compared to favipiravir, lopinavir, and hydroxychloroquine (−10.76, −10.72 and −10.35 kcal/mol, respectively) as shown in Table 2. Luteolin interacts with ASN 544, and GLY 545 amino acid residuals of 6VYB through hydrogen bonds and the distance of these bonds are 2.32, and 2.16 Å, respectively. It also interacts with ALA 522, and HIS 519 amino acid residues through Van der Walls interaction. The ASN 969, and ASP 867amino acid residuals of 6VYB interacts with mundulinol through hydrogen bonds with a distance of 2.33, and 2.39 Å, respectively. It also interacts with PHE 970 amino acid residues through Van der Walls interactions. Basically, in the state of 6W63, the binding affinity of the studied compounds is comparatively higher than those of 6VYB. When we compare the docking results of docked ligands with 6W63 and 6VYB proteins, we have noticed that luteolin and mundulinol compounds have higher binding energies than currently used antiviral drugs in both cases (6W63 and 6VYB). Whilst derrisin failed to have higher binding energy than those of favipiravir, lopinavir, and hydroxychloroquine, derrisin is excluded from the promising drugs list as antiviral against SARS CoV-2. Based on this, only luteolin and mundulinol compounds exhibited the ability to work as dual inhibitors against 6VYB and 6W63 because both had binding affinities higher than those of selected drugs favipiravir, lopinavir, and hydroxychloroquine (as anti SARS-Cov-2 drugs) in two cases (6VYB and 6W63).

ADMET Properties and Drug-Likeness Model Score
All the values obtained from the ADMET analysis for our studied compounds are presented in Table 3. Lipinski's rule of five is a widely used method for drug development [49]. All compounds in our study complied with Lipinski's rules of five. Based on BBB permeability, the compounds were non-BBB permeable. Additionally, all the studied compounds have high GI absorption, especially apigenin, luteolin, and mundulinol compounds with optimum solubility. The solubility and Caco-2 permeability are used to evaluate the drug absorption. The solubility ranges from −5.86 to −1.49 indicating the moderate to the high solubility of all studied compounds. The high-water solubility can be due to different hydroxy groups present in each compound capable of forming bonds with water molecules. The human colonic adenocarcinoma cell line Caco-2 permeability of all studied compounds have high permeation. The gastrointestinal absorption of all studied compounds. All ADMET properties of all studied compounds are in a good agreement of those calculated by Vijayakumar et al. [50]. Vijayakumar [51] studied the activity of some flavonoids compounds against chikungunya virus replication using molecular docking and DFT-based approach. They calculate the Pharmacokinetics and pharmacological properties of some flavonoids compounds, they elucidate that Apigenin and Luteolin have high solubility and GI absorption. Additionally, the toxicity and carcinogenicity of both compounds were also predicted for in order to determine the safe use of the selected compounds during in vivo analysis. Apigenin and Luteolin showed non-toxic and non-carcinogenic nature. Rasool et al. [52] studied the inhibitory potential of selective phytochemicals against Mpro of 2019-nCoV by computer-aided approaches, they study ADMET properties of Mundulinol, Derrisin, Luteolin, and Apigenin. They found that Mundulinol, Luteolin, Derrisin, and Apigenin have high soulibilty and GI absorption. Additionally, four compounds exhibited the properties of non-BBB permeability. Finally, four compounds exhibited non-toxic and non-carcinogenic nature. In a collective sense from all cited published paper in this section that, the ADMET profile of Luteolin, Mundulinol, and Apigenin compounds is satisfactory and hence suitable for in silico studies with SARS-CoV-2 proteins. Finally, estimations of the toxicity and carcinogenic properties, suggesting that all studied compounds are non-toxic and non-carcinogenic. So, all compounds have suitable ADMET profiles, especially compounds such as apigenin, luteolin, and mundulinol. The druglikeness model score for the selected compounds was predicted from the Molinspiration server (http://www.molinspiration.com accessed on 10 June 2021) and the output values are shown in Table 4. We excluded the compounds with zero or negative value of druglikeness score from nominated compounds to be drugs (derrisin and rhamnetin), while the highest drug-likeness score was found to be 0.98 for apigenin followed by luteolin (0.95), mundulinol (0.92), and finally naringenin (0.90). Joining together the docking results, ADMET estimations, drug-likeness results suggest that compounds, luteolin, and mundulinol are promising to be developed as leads drugs with good pharmacokinetic properties. Therefore, we nominated these two compounds for further investigations via running MD simulations.

Molecular Dynamic Simulation (MD)
MD simulation mimics the structural perturbations and actual movement of a receptor molecule within corresponding biological conditions. Our docking poses are validated through the RMSD, a well-established parameter for protein stability and equilibration, and examine the average behavioural pattern for the entire protein structure during the MD simulation analysis.

Root Mean Square Deviation (RMSD)
The RMSD score determines the direct changes in a protein from its primary coordinates. The RMSD values for the protein backbone coupled with its potential inhibitors were estimated using the initial structure as a frame reference (0 to 100 ns). RMSDs of free 6W63 receptor, mundulinol + 6W63, and luteolin + 6W63 complexes are shown in Figure 5. Additionally, RMSDs of free 6VYB receptor, mundulinol + 6VYB, and luteolin + 6VYB complexes are presented in Figure 6. The RMSD values gradually augmented from 0 to 20 ns and retained the values constant after reaching an equilibration throughout the simulation period for all the studied systems. The RMSD values for all studied complexes were fluctuating between 20 to 40 ns indicating that the compounds were familiarizing a new conformation within the binding pocket. The average RMSD values for the last 20 ns for the systems, were 0.430 ± 0.06, 0.435 ± 0.09, 0.440 ± 0.14, 0.525 ± 0.03, 0.530 ± 0.07, and 0.540 ± 0.09 for 6W63 receptor, mundulinol + 6W63, luteolin + 6W63, 6VYB receptor, mundulinol + 6VYB, and luteolin + 6VYB complexes, respectively. Lower RMSD values of all studied complexes indicate that mundulinol and luteolin inhibitors are stable inside 6W63 and 6VYB and provide a significant base for our study.

Root Mean Square Fluctuation (RMSF)
RMSF was calculated in respect of the backbone atom per amino acid residues and the fluctuations at residue level were demonstrated by the RMSF plot. Figure 7 shows RMSFs of 6W63 receptor, mundulinol + 6W63, and luteolin + 6W63 complexes. Figure 8 shows RMSFs of 6VYB receptor, mundulinol + 6VYB, and luteolin + 6VYB complexes. Both the free receptors (6W63 and 6VYB) and the complexes with low RMSF values showed the nearly same behaviour of residue fluctuation report. So, the RMSF plot confirmed that the binding of all studied compounds is stable inside the receptors also, there is no major consequence found during the simulation analysis regarding the flexibility of the receptor. To gain better understandings of the receptor flexibility, the time average of RMSF values of the 309 amino acids of 6W63 and 975 amino acids of 6VYB receptors in the absence and presence of the studied compounds during the simulation period was computed. In the state of 6W63, the RMSF values for the complexes reported less fluctuation for the catalytic dyad residues (His41 and Cys148). The average RMSF values were 0.33 ± 0.04, 0.40 ± 0.07, and 0.47 ± 0.07 Å for 6W63 receptor, mundulinol + 6W63, and luteolin + 6W63 complexes, respectively. In the state of 6VYB, the RMSF values suggested that the following residues (HIE-41, LEU-141, and GLN-189) showed less fluctuation in the system. The average RMSF values were 0.50 ± 0.06, 0.53 ± 0.04, and 0.55 ± 0.05 Å for 6VYB receptor, Mundulinol + 6VYB, and Luteolin + 6VYB complexes, respectively.

Hydrogen Bonds Analysis
Hydrogen bonds analysis is used to gain more insights into the molecular recognition, molecular interactions, and selectivity of the studied compounds inside the receptors. During MD simulation, the protein-ligand interactions obtained from the changes in the secondary structures were controlled by the analysis of hydrogen bonds. There are different conformations obtained from MD simulation in which each conformation of a receptor has its specific model of interaction with the ligand. We analysed the number of hydrogen bonds formed of MD simulations of all trajectories for the complexes. As showed in Figure 9, the number of hydrogen bonds formed during MD simulation time for mundulinol + 6W63 and luteolin + 6W63 complexes, the average number of conformations were nearly three hydrogen bonds for mundulinol − 6W63 complex and two hydrogen bonds for luteolin − 6W63 complex. Figure 10 represents the mean number of hydrogen bonds generated during MD simulations time for mundulinol − 6VYB complex and luteolin − 6VYB complexes was higher than the number of established hydrogen bonds with SARS-CoV-2 main protease. Further, the number of hydrogen bonds of mundulinol − 6VYB is six while the number of hydrogen bonds of luteolin − 6VYB is five. These results showed that the mundulinol and luteolin compounds possessed a large number of hydrogen bonds with 6VYB during the simulation than 6W63. So, mundulinol and luteolin were capable of maintaining a strong interaction with the binding pocket of 6VYB and 6W63 through the simulation period.

MM-PBSA Binding Free Energy
MM-PBSA calculations for two selected complexes were estimated to analyse the biophysical basis [53] of the molecular recognition of Mundulinol and Luteolin as dual targets inhibitors. It is evident from Table 5 that the electrostatic interactions (∆E Electrostatic ) values of Luteolin are more favoured than Mundulinol. Further, the available data in Table 5 reveals that ∆E Van derWaal values (intermolecular van der Waals) are similar for Luteolin and Mundulinol when they are targeted SARS-CoV-2 main protease M pro (PDB ID: 6W63) and Mundulinol is a little more favoured than Luteolin when they bind to SARS-CoV-2 spike protein Sp (PDB ID: 6VYB). The Lower binding free energy determines the better receptor and ligand interaction. Based on mentioned parameters, Table 5, Luteolin compound showed better binding free energy when it binds to 6VYB and 6W63 when compare with Mundulinol. Additionally, the binding free energies of both compounds in the state of 6VYB are greater than those of 6W63. This indicated the compounds as stable and had a negligible effect on the flexibility of the 6VYB and 6W63 during the simulation periods and reinforced the significance of our study. ∆E binding can offer a more realist view concerning the affinities of Mundulinol and Luteolin against the selected targets.

Ligand and Protein Preparation
The structures of studied flavonoid derivatives were downloaded from the Pub-Chem database and saved in SDF format. The energy minimization and structural optimization of the compounds were performed by using the Molecular Mechanics 2 (MM2) force field. Next, the Autodock tool was applied for saving the energy minimized structures in pdbqt format. The 3D crystal structure of spike protein (PDB ID: 6VYB) and main protease (PDB ID: 6W63) of SARS CoV-2 were obtained from the protein data bank (https://www.rcsb.org/ (accessed on 10 June 2021).) [54]. Elimination of water and the small molecules from crystal structures was accomplished via BIOVA Discovery Studio software [52]. Following, the joining of the polar hydrogen atoms and processing the atomic charges were completed by Kollman and Gasteiger method [55]. Finally, the 3D structure of proteins was kept in pdbqt format.

Molecular Docking Analysis and Binding Energy Estimation
The inhibitory role of our studied flavonoid compounds against the spike protein (PDB ID: 6VYB) and main protease (PDB ID: 6W63) of SARS-CoV-2 was analysed by the molecular docking study. The molecular docking study, in respect of binding of these flavonoids with receptors and their reactivity, following methods applied in various studies [58][59][60]. All the parameters were set at default values for docking analysis [53]. The docking grid covered all the active sites of the spike protein and main protease receptors with a grid size of 70 Å × 70 Å × 70 Å and a spacing value of 0.375 Å. The docking analysis used the Lamarckian Genetic Algorithm (LGA); all the calculations for docking analysis were carried out using the autodock 4.2 software [61]. Further evaluation of the docking output was completed using the discovery studio program. After the docking is completed, the best pose was chosen based on the lowest binding energy. Later on, we aligned our docked poses with that of co-crystallized structure, and finally RMSD lower than 1.0 Å was calculated.

Molecular Dynamic (MD) Simulation
Molecular dynamic simulation is broadly utilized to evaluate the biological molecules at the atomic level [62,63] including structural stability and conformational changes. The compounds luteolin and mundulinol exhibit the lowest docking energy with the main protease and spike protein of SARS-CoV-2. The lowest energy of compounds denotes favourable binding affinity [64], therefore it can be said that these the compound could be treated as an inhibitor. Herein, top-ranked compounds bound complexes obtained from molecular docking analysis (mundulinol + 6W63, luteolin + 6W63, mundulinol + 6VYB, and luteolin + 6VYB complexes) were subjected to NAMD software for MD simulation [65]. The temperature and pressure were set to 303 K and 1 atm, respectively. All complexes were solvated using the transferable intermolecular potential water molecules (TIP3P) model. The CHARMM force field was used for a short energy minimization (1000 numsteps) by the Steepest Descent method. Following, the MD simulations were accomplished in two steps: equilibration under NVT ensemble for 5 ns timescale, and production for 100 ns timescale under the NPT ensemble. All theoretical background of MD simulations was performed according to Al-Anazi et al. [66].

MM-PBSA Calculation
The free energy for binding (∆G bin ) was estimated via the calculation of the Free Energy (CaFE) program exploiting the Molecular Mechanics-Poisson Boltzmann Surface Area (MM-PBSA) method [67,68]. It handles simulation trajectories produced from numerous force fields and is enforced by VMD and NAMD. Three distinct energy terms are estimated through the MM-PBSA method that includes energy variations between the complexes and keeps the receptor and ligand molecule apart from each other in the gas phase from NAMD. Then, the APBS program was used to measure the polar solvation free energy [69], Lastly, the difference of solvent-accessible surface area (SASA) and the free energy of nonpolar solvation was also estimated. The energies needed for binding of the ligand-receptor complexes were calculated throughout the whole MD simulation (100 ns).

Conclusions
Current research endeavours against COVID-19 involve numerous antiviral and conventional drugs to fight this pandemic, yet no ultimate therapeutic interventions are available so far. Using the molecular docking analysis, we found the interaction of luteolin and mundulinol with SARS-CoV-2 main protease, spike protein has the highest binding energy when compared to all studied compounds and some proposed antiviral drugs. So, the molecular docking study revealed that luteolin and mundulinol compounds can act as potential inhibitors for Mpro and Sp. The ADME/Tox properties strongly suggest all the studied compounds were non-toxic and non-carcinogenic. Finally, the RMSD value of the receptor-ligands complexes calculated for both receptors exhibited stability at around 0.6 Å and RMSD of two compounds (luteolin and mundulinol) complexed inside 6W63 and 6VYB remained at a favourable range of 0.3 Å and retained its stability throughout the simulation process. The free receptor and backbone atoms of the complex exhibit comparable RMSF values, showing the stability of the studied compounds inside 6W63 and 6VYB receptors. On the whole, the favourable effects of luteolin and mundulinol led to a conclusion that these may be considered as promising treatment strategies against COVID-19 after further clinical trials.