A threefold approach including quantum chemical, molecular docking and molecular dynamic studies to explore the natural compounds from Centaurea jacea as the potential inhibitors for COVID-19

Abstract In the current report, we studied the possible inhibitors of COVID-19 from bioactive constituents of Centaurea jacea using a threefold approach consisting of quantum chemical, molecular docking and molecular dynamic techniques. Centaurea jacea is a perennial herb often used in folk medicines of dermatological complaints and fever. Moreover, anticancer, antioxidant, antibacterial and antiviral properties of its bioactive compounds are also reported. The Mpro (Main proteases) was docked with different compounds of Centaurea jacea through molecular docking. All the studied compounds including apigenin, axillarin, Centaureidin, Cirsiliol, Eupatorin and Isokaempferide, show suitable binding affinities to the binding site of SARS-CoV-2 main protease with their binding energies -6.7 kcal/mol, -7.4 kcal/mol, -7.0 kcal/mol, -5.8 kcal/mol, -6.2 kcal/mol and -6.8 kcal/mol, respectively. Among all studied compounds, axillarin was found to have maximum inhibitor efficiency followed by Centaureidin, Isokaempferide, Apigenin, Eupatorin and Cirsiliol. Our results suggested that axillarin binds with the most crucial catalytic residues CYS145 and HIS41 of the Mpro, moreover axillarin shows 5 hydrogen bond interactions and 5 hydrophobic interactions with various residues of Mpro. Furthermore, the molecular dynamic calculations over 60 ns (6×106 femtosecond) time scale also shown significant insights into the binding effects of axillarin with Mpro of SARS-CoV-2 by imitating protein like aqueous environment. From molecular dynamic calculations, the RMSD and RMSF computations indicate the stability and dynamics of the best docked complex in aqueous environment. The ADME properties and toxicity prediction analysis of axillarin also recommended it as safe drug candidate. Further, in vivo and in vitro investigations are essential to ensure the anti SARS-CoV-2 activity of all bioactive compounds particularly axillarin to encourage preventive use of Centaurea jacea against COVID-19 infections.

SARS-CoV-2 Xu et al., 2020). The M pro is a key enzyme of SARS CoV-2 with molecular mass of 33796.8 Da and plays a pivotal role in processing of two polyproteins pp1a and pp1ab into 16 NSPS (non-structural proteins). These NSPS further involve in the synthesis of subgenomic RNAs, these RNAs encode for accessory proteins and four different structural proteins (envelope (E), membrane (M), spike (S), and nucleocapsid (N) proteins). Therefore M pro , plays a vital role in the proteolytic maturation and life cycle of SARS CoV-2 (Ramajayam et al., 2011, Ren et al., 2013. Thus, inhibition of M pro can be a potential target for the inhibition of SARS CoV-2 life cycle at transcription or translation phase. Moreover, there is no homolog of M pro in human genome thus inhibitors of M pro will be safe for human (Hayden et al., 2003, Yang et al., 2005 which makes it a supreme antiviral drug target . Now a day computational biology is playing a key role in drug discovery. Molecular docking has been used to check the interactions of ligand (drug) with the binding sites of targeted protein (receptor) (McConkey et al., 2002). It deals with all possible factors involved in drug discovery such as identification of hit compound, optimization of most appropriate molecule and virtual screening (Jorgensen, 2004, Kitchen et al., 2004. Moreover, in silico studies has proved to be very useful in drug discovery in saving resources in terms of time as well as money (Murgueitio et al., 2012). Besides this, in silico studies increases the potential to discover new therapeutic drugs while reducing extensive lab work. Nonetheless, these studies use with care and reported drug targets should test in lab experiments to get benefit of in silico studies.
Designing drugs with maximum selectivity to the target molecule is the main idea of drug discovery (Gertsch, 2011). Medicinal plants are well known among the key sources to provide new pharmaceuticals with maximum selectivity (Balunas and Kinghorn, 2005;Chiarello, 1995). Herbal medicines have a number of advantages over modern synthetic drugs, including fewer side effects, affordability, maximum selectivity with the receptor and cost effectiveness encouraging the discovery of plant based drugs (Chakraborty, 2018;Sharifi-Rad et al., 2018). Plants produce different types of secondary metabolites such as phenols, flavonoids, and alkaloids. Plants may synthesize these metabolites for growth or therapeutic purposes Chakraborty, 2018, Kumar et al., 2020. Moreover, different plants such as Centaurea aegyptiaca, Centaurea alba etc. have already been investigated due to their antiviral effects (Bakr et al., 2016;Politeo et al., 2019). Centaurea jacea commonly known as brown knapweed is a perennial herb of family Asteraceae. Sporadically it has

Introduction
Viral infections are continued to emerge and symbolize a major problem to public health. In last 20 years, numerous viral epidemics such as the SARS-CoV (severe acute respiratory syndrome coronavirus) in 2002 to 2003, H1N1 pandemic in 2009 and MERS-CoV (Middle East respiratory syndrome coronavirus) was first discovered in Saudi Arabia in 2012, have been reported (Perlman and Netland, 2009;Chan et al., 2013). Previous studies have described that CoVs (coronaviruses) are an etiologic agent of various infections especially respiratory and digestive disorders in mammals, reptiles, avian species and humans (Malik et al., 2020). In the last weeks of 2019, a viral infection COVID-19 caused by a novel virus named as SARS-CoV-2 (severe acute respiratory syndrome-coronavirus 2), a new strain of CoVs has been reported in Wuhan, China and spread rapidly throughout the whole world . SARS-Cov-2 has been identified as an enveloped, non-segmented RNA virus with 29.5kb genome (Adem et al., 2020). Genome sequencing of SARS-CoV-2 showed that SARS-CoV-2 is 96.2% identical to a bat CoV and the RNA genome of SARS-CoV-2 is 79.5% identical with the RNA genome of SARS-CoV Wu et al., 2020;Zhou et al., 2020).
A recent report of WHO described that epidemic COVID 19 increases rapidly to more than 216 countries, till 03 April 2021, worldwide 129,215,179 cases and 2,820,098 deaths were reported due to COVID-19 (WHO, 2021). Fast spreading rate and no specific antiviral drug for COVID-19 are few major reasons for such outspread of this pandemic. Currently no specific antiviral drugs or vaccines are available for the remedy of COVID-19. However, recent research suggested that cotreatment of COVID-19 patients with certain previously used antiviral drugs can yield optimistic results (Chang et al., 2020). In emergency situations combination of Lopinavir/ Ritonavir (HIV drugs) and α-interferon has been used for the treatment of COVID-19 patients, but therapeutic effects remain very restricted and there can be side effects too (Cao et al., 2020). Other described antiviral compounds for human CoVs include neuraminidase inhibitors, nucleoside analogues, remdesivir, tenofovir disoproxil, umifenovir (arbidol) and lamivudine . In this need of hour specific targeted anti-SARS-CoV-2 drugs with safety and efficacy are urgently needed.
Palavras-chave: Covid-19, Apigenina, Centaurea jacea, química quântica, docking molecular, dinâmica molecular. been used in different folk medicines for the treatment of dermatological complaints and fever. Its chloroform extract contains different secondary metabolites such as centaureidin, isokaempferide, axillarin, eupatorin, hispidulin, apigenin and cirsiliol (Forgo et al., 2012). Secondary metabolites of Centaurea jacea are bioactive molecules and used in folk medicines for centuries due to their anti-inflammatory, antioxidant, antiviral and antibacterial activities (Yan et al., 2017). Thus, in the present study, we have studied different bioactive molecules of Centaurea jacea against M pro by in silico molecular docking study. We plan not only to disclose the most effective bioactive compound of Centaurea jacea as a potent inhibitor of M pro but also explore several structure-activity relations by studying intermolecular interactions between potent inhibitors and of M pro protein.

Computational Methods
AutoDock Vina (ADV) was used to perform all the docking analysis because (a) ADV offers more accuracy in analyzing protein-ligand interaction as compared to AutoDock 4.2 (b) it takes shorter running time due to its several core processors (c) ADV provides more accuracy for ligand analyzing more than twenty rotatable bonds (Narkhede et al., 2020). Moreover, autodock MGL tools 1.5.6 and Pymol were used to obtain PDBQT format of ligand and protein. Finally, docking results were visualized by using free source BDSV (BIOVIA Discovery Studio Visualizer (Narkhede et al., 2020). Docking was done as a blind docking (blind docking refers to the use of a grid box (40 number of points in X, Y and Z dimensions with -26.283 12.599 58.966 X, Y and Z center respectively and spacing 0.375 Å) which is large enough to encompass any possible ligand-receptor complex) using ADV.

Receptor preparation
The M pro of SARS-CoV 2 (PDB ID: 6LU7)  was used as the rigid receptor. AutoDock MGL Tools 1.5.6 was used to prepare the receptor protein. Preparation of receptor involves removal of previously attached ligand with M pro , removal of all water molecules, addition of polar hydrogen atoms, addition of kollman charges and finally to save it in PDBQT format.

Ligand preparation
The 3D sdf formats of all ligands were downloaded from pubchem (Kim et al., 2021). PDBQT format of all ligands were obtained by using Auto Dock MGL Tools 1.5.6. The docking was performed using exhaustiveness value of 8, remaining parameters of the software were sustained as a default and all bonds of the ligand were allowed to move freely, taking receptor as a rigid. The final docking results were visualized through Discovery Studio Visualizer 2.5.

ADME Analysis
ADME (Adsorption, Distribution, Metabolism and Excretion) is important to investigate the pharmacodynamics of proposed compound. A website SWISS-ADME (Daina et al., 2017) allows the user to enter SMILES or draw structure of proposed drug or ligand from Pubchem and offers the parameters such as lipophilicity (iLOGP, XLOGP3, WLOGP, MLOGP, SILICOS-IT, Log P0/w), drug likeness rules (Lipinski, Ghose, Veber, Egan and Muegge), water solubility-Log S (ESOL, Ali, SILICOS-IT) and Medicinal Chemistry (PAINS, Brenk, Leadlikeness, Synthetic accessibility) methods are analyzed (Daina et al., 2017). SMILES of all the ligands were obtained from Pubchem and entered into the search bar of SWISS-ADME and analyzed.

Toxicity prediction
Prediction of toxicology is important to predict the quantity of tolerability of molecule before being used by human or animal model. The pkCSM is an online website in which ligand can be drawn virtually or can be studied by submitting SMILES of the ligand. The website provides the details of toxicology in multiple fields such as maximum tolerable dose by human, AMES Toxicity, hERG-I and hERG-II inhibitor, LOAEL, LD50, Skin Toxicity, Hepatotoxicity, T. pyriformis and Minnow toxicity (Pires et al., 2015). SMILES of ligands were searched from PubChem and submitted into online database pkCSM and toxicity mood was selected.

Structural chemistry of ligands
It pertinent to discuss the ligands chemistry before understanding docking results. The geometrical structures of all six ligands were optimized using density functional theory methods with B3LYP/6-311G* level of theory as implemented in Gaussian 16 suit of programs (Frisch et al., 2009). All the optimized ligands showed the planner lowest energy geometries with slightly out of plane methoxy groups. The molecular electrostatic potential (MEPs) diagrams are calculated to see the reactivities of the ligands. The MEP provides many comprehensive intuitives about the distribution of electrostatic charges on the total density surface of optimized ground state geometries of ligands as shown in Figure 1. For instance, a maximum negative region is preferred site for electrophilic attack as it is indicated with red surface. In our studied ligands, the hydroxyl and carbonyl groups are showing negative potential regions. On the other hands, the maximum positive region that is preferred site for nucleophilic attack indicated as blue color on MEP surface which are mostly acidic and/or polar hydrogen atoms of ligands. Furthermore, we have also illustrated the HOMO and LUMO orbitals of all ligands. The frontier molecular orbitals especially HOMO and LUMO play a significant role to determine the reactivity of a chemical compound. The HOMO and LUMO orbitals are shown in Figure 1 with their respective orbital energies. A careful analysis of HOMO orbitals shows that ligands Axillarin, Apigenin, Isokaempferide and Centaureidin possess somewhat uniform distributions of their valence electronic clouds as compared with other counterparts. Furthermore, we have calculated the HOMO-LUMO energies and their energy gaps (∆ HL ) which are very crucial to define their reactivities. A graphical illustration of ∆ HL values were given in Figure 2, which showed that axillarin ligand possessed the lowest ∆ HL value as compared to the other studied ligands. Usually, a lower ∆ HL value shows higher kinetic reactivities of ligands (Aihara, 2000). Additionally, it is important to note that the number of hydroxyl, carbonyl and methoxy groups is also very vital for establishing durable interactions with receptor proteins.

Binding energy and inhibition constant
Docking analysis was conducted triplicate to get a clear demonstration of the mode of action of various bioactive constituents of Centaurea jacea against M pro . For each constituent of Centaurea jacea nine interactions were generated and the one was selected with lowest binding energy. Our findings depicted that all the constituents of Centaurea jacea can be a potential inhibitor of M pro .
The binding pocket of M pro comprises of two chains, Chain A and C which may play a pivotal role in interaction with ligand. Ligand may interact either with chain A or chain C depending upon the obtainability of atoms for specific interactions (Table 1). Table 1 shows the binding energies, inhibition constant and interacting amino acid residues obtained from the docking of M pro against different bioactive components (Apigenin, Axillarin, Centaureidin, Cirsiliol, Eupatorin, Isokaempferide) of Centaurea jacea. Apigenin, Axillarin, Centaureidin, Cirsiliol, Eupatorin and Isokaempferide were found as a best potential inhibitor of M pro , all of these ligands were bind with the chain A of M pro . The binding energies acquired from the docking of M pro with Apigenin, Axillarin, Centaureidin, Cirsiliol, Eupatorin and Isokaempferide were -6.7 kcal/mol, -7.4 kcal/mol, -7.0 kcal/mol, -5.8 kcal/mol, -6.2 kcal/mol and -6.8 kcal/mol respectively (see Table 1 and Figure 3).
Inhibition constant obtained from the binding energies of M pro in complex with Apigenin, Axillarin, Centaureidin, Cirsiliol, Eupatorin and Isokaempferide were 11.77 µM, 3.59 µM, 7.08 µM, 54.09 µM, 27.46 µM and 9.94 µM respectively see table 1. Figure 4 shows the direct relation between the binding energy and inhibition constant of all the ligands, it revealed that axillarin has lowest inhibition constant 3.59 µM with lowest binding energy -7.4 kcal/mol similarly Cirsiliol has the highest inhibition constant 54.09 µM with the highest binding energy -5.8 kcal/mol.

Protein-ligand interactions
Discovery studio visualizer was used to visualize the interactions between different constituents of Centaurea jacea and M pro . Binding of ligands in different binding sites of M pro are shown in Figure 5, protein is presented as a surface while ligand is displayed as stick model, moreover green and pink surface of the protein show acceptor and donor regions of the protein respectively Figure 5. The amino acids involve in ligand protein interactions are displayed with ligands as yellow stick model with different amino acids surrounding them. Green dashes show hydrogen bond between protein and ligand. Different amino acids involve in protein ligand interactions are displayed as stick of different colors and labeled by red color as in Figure 6. Figure 7 shows a 2D diagram of protein ligand Results of our study revealed that apigenin forms 3 hydrogen bonds with LYS137 (2.32 Å), LYS137 (2.40 Å) and GLU290 (2.82 Å) protein residues, 1 hydrophobic interaction with LEU286 protein residue and 2 electrostatic interactions with LEU286 amino acid residue of M pro . Axillarin form 5 hydrogen bonds with HIS41 (2.91 Å), CYS145 (3.56 Å), GLU166 (2.82 Å), HIS172 (3.02 Å) and PHE140 (2.22 Å) protein residues and 5 hydrophobic interactions with MET49, HIS41, LEU141, MET49, and MET165 amino acid residues of M pro . The most striking property of axillarin is its interactions with CYS145 and HIS41 which are the most crucial binding sites of M pro (Das et al., 2021). Centaureidin form 6 hydrogen bond with GLY143 (1.95 Å), GLY143 (2.54 Å), GLU166 (2.42 Å), ARG188 (2.53 Å), LEU141 (2.86 Å) and THR26 (3.59 Å) amino acid residues, 4 hydrophobic interactions with LEU27, CYS145, CYS145 and MET165 while one electrostatic interaction with HIS41 amino acid residues of M pro . Cirsiliol potentially form 7 hydrogen bonds with LYS102 (2.48 Å), LYS102 (2.24 Å), LYS102 (2.87 Å), ASN151 (2.82 Å), ARG105 (2.25 Å), ASP153 (2.10 Å) and SER158 (3.24 Å) protein residues while 3 hydrophobic interactions with VAL104 protein residue of M pro as in Table 1      H-bonds play a vital role in stabilizing the structures of all proteins due to their involvement in the secondary structure elements such as beta sheets and alpha helix. Moreover H-bond interactions are important elements for ligand binding affinity. Specificity of ligand binding is mainly determined by the number of intermolecular H-bonds involving crucial residues present as an acceptor and donor in protein and ligand structures. These H-bond networks played a significant role in supporting the binding affinity between ligand and protein (Bitencourt-Ferreira et al., 2019). As Figure 8 shows that in our study most of the interactions between receptor-ligand complexes are hydrogen bond interactions so we also plot a graph between the H-bond distances on Y-axis while protein residues according to their respective ligand was plotted on X-axis Figure 8. The shorter the bond distance will result in stronger bond between receptor-ligand complexes. Figure 9 shows that in our study H-bond between centaureidin and GLY143 residue of M pro have minimum bond distance 1.98 Å while H-bond between eupatorin and TYR239 residue of M pro have maximum bond distance 3.76 Å. H-bond distance between axillarin and respective M pro residues ranges in 2.22-3.56 Å. Table 1 shows the bond distance of all protein residues interacting with different ligands.

Root Mean Square Deviation (RMSD) Study
Among the studied complexes, we choose axillarin-M pro complex for further molecular dynamic studies owing to its highest binding affinity and good intermolecular interactions. The NAMD (Phillips et al., 2020) with CHARMM force-field (Best et al., 2012) was used to performed minimization followed by dynamic run while trajectory, RMSD, RMSF and other analysis were done by VMD program (Humphrey et al., 1996) as developed at NIH Center for Macromolecular Modeling and Bioinformatics, Theoretical and Computational Biophysics group. The MD simulations of axillarin-M pro complex and M pro protein were performed for 60 ns (6×10 6 femtosecond) time scale. Further details of molecular dynamic calculations can be found in the supporting information of the article. The Figure 10 shows the RMSDs of axillarin-M pro complex and M pro protein in overlapping fashion. A careful analysis of RMSDs indicates that the RMSD of M pro protein (blue curve) shows some fluctuations (between 0 to 2.7 Å) for the first 0 -7 ns and after that its RMSD becomes relatively consistent which only fluctuates between 1.5 to 2.5 Å throughout the remaining time scale. While on the other hands, the RMSD of axillarin-M pro complex shows a different behavior where at the start, it shows significant variations during 1 st nanosecond of trajectory and after that its variations remain between 1.2 to 2.7 Å till 40 ns. From 40 to 50 ns, the RMSD of axillarin-M pro complex is upset more, which finally show some convergence from 50 to 60 ns by constraining the fluctuations between 1.8 to 2.4 Å. The RMSDs analysis indicates that the docking of ligand has disturbed the RMSD of M pro protein (blue curve) especially at the start and from 40 to 50 ns and then show reasonable stability from 50-60 ns.

Root Mean Square Fluctuations (RMSF) of M pro residues
The RMSF is a useful analysis to describe the local changes along the chain in protein residues. We have   Synthetic accessibility for axallrin is found to be 3.46. ADME properties for rest of the compounds are given in Table 2.
Molar refractivity of the drug compounds endorses that the compound is permeable for specific membranes and can persist even in the midst of weak or strong solventsolvent and solute-solvent interactions. Lipophilicity of a compound tells us about the oral, intestinal and sublingual absorption of the compound. Water solubility values of a compound shows that the compound is freely soluble in water or not (Enmozhi et al., 2021).

Toxicity prediction
Toxicity prediction analysis for different constituents of Centaurea jacea was carried out by using an online website, pKCSM (http://biosig.unimelb.edu.au/pkcsm/ prediction). Results of the analysis shown in Table 3, all of the components of Centaurea jacea don't show AMES toxicity, don't have ability to induce hepatotoxicity and do not cause skin sensitivity. Maximum tolerable dose for human of our best docked compound axillarin is about 0.506 log mg/kg/day, acute oral rat toxicity (LD50) and chronic oral rat toxicity is 2.35 mol/kg and 2.50 log mg/ kg_bw/day respectively. Moreover, 0.30 log mg/L and 2.979 log mM of axillarin can cause T. pyriformis and Minnow Toxicity, respectively. Molecular docking, ADME analysis and toxicity prediction results of axillarin and M pro complex are appreciable in silico, therefore, in vitro, in vivo and clinical studies should be done for the development of novel antiviral drug against novel SARS CoV-2.
Aforementioned results of molecular docking show the axillarin as an ideal inhibitor of M pro especially interaction calculated and plotted the RMSFs of M pro protein and RMSF of M pro while it is docked with axillarin to see the effect of docking as shown in Figure 11. In the RMSF graphs, the peak areas indicate the residues which fluctuates the most during the simulation. It is commonly observed that tail areas with C-and N-terminals show larger fluctuations as compared to the other parts. The overlapping of RMSFs of axillarin-M pro complex and M pro protein shows that the initial three residues show maximum fluctuations both in M pro protein and its complex with axillarin, which is owing to their primary structural configurations in M pro protein. A careful analysis of both RMSFs illustrates that there is reasonable stability after docking of axillarin with M pro protein because overall no abnormal and/or significant fluctuations are seen in RMSFs of M pro protein and axillarin-M pro complex (see the red curve).

Predicted ADME (absorption, distribution, metabolism and excretion) properties of ligands
To predict the drug likeness, the ADME properties for various constituents of Centaurea jacea were calculated (Table 2). Molar refractivity of our best docked compound axallrin is 89. The LogP values of ADME analysis for axillarin represented its lipophilic nature with consensuses LogP values 2.30. Water solubility value of ADME analysis for axillarin is -3.81 which suggested that axalllrin is moderately soluble in water. Axillarin is considered as a good candidate of drugs for oral administration because it is highly absorbable in GIT (gastrointestinal tract). Axillarin clearly fulfilled the Lipinski's rule of five (Lipinski, 2004). potential of different compounds as potential inhibitor for COVID-19 M pro from Centaurea jacea. All the studied compounds such as apigenin, axillarin, Centaureidin, Cirsiliol, Eupatorin and Isokaempferide, could strongly bind to the binding site of SARS-CoV-2 main protease with binding energies -6.7 kcal/mol, -7.4 kcal/mol, -7.0 kcal/mol, -5.8 kcal/mol, -6.2 kcal/mol and -6.8 kcal/mol respectively. Inhibition constant obtained from the binding energies of M pro in complex with Apigenin, Axillarin, Centaureidin, Cirsiliol, Eupatorin and Isokaempferide were 11.77 µM, 3.59 µM, 7.08 µM, 54.09 µM, 27.46 µM and 9.94 µM respectively. Our results proved that studied bioactive constituents of Centaurea jacea strongly bind with the of axillarin with the most important catalytic residues of M pro proves it as a potential inhibitor of M pro . Moreover, druglikness and toxicity prediction analysis also revealed the use of axillarin as a safe and specific antiviral drug against SARS-CoV-2. As it's an in-silico study, therefore in vivo and in vitro studies should be done in future for the development of axillarin as an anti-SARS-CoV-2 drug.

Conclusions
Thus, we have successfully performed molecular docking and molecular dynamic studies to explore the  different binding sites of M pro of SARS-CoV-2, which is essential for proteolytic maturation of SARS-CoV-2. Among all the docked compounds, the axillarin was found to be the most effective inhibitor of M pro based on its lowest binding affinity of -7.4 kcal\mol, inhibition constant 3.59 µM and its maximum hydrophobic, electrostatic and hydrogen bond interactions. Furthermore, the molecular dynamic calculations also showed significant insights into the binding effects of axillarin with M pro of SARS-CoV-2. The RMSD and RMSF calculation indicates the stability and dynamics of complex in aqueous environment. Besides this, the results of quantum chemical study for the ligands also show that axillarin is most reactive among all docked compounds due to lowest ∆ HL . In addition to molecular docking study druglikness prediction shows that axillarin clearly follow the Lipinski rule of drug likeness, moreover toxicity prediction also proves it as a safe drug for human. Results of Molecular docking, druglikeness and toxicity prediction for all the studied compounds of C. jacea particularly for axillarin encourage further in vivo, in vitro and clinical studies to prove its inhibitory effect against M pro of SARS-CoV-2 of axillarin.