In silico screening of potential inhibitors from Cordyceps species against SARS-CoV-2 main protease

Abstract Coronavirus disease 2019 (COVID-19) is a result of a retroviral infection of SARS-CoV-2. Due to its virulence and high infection rate, it is a matter of serious concern and a global health emergency. Currently available COVID-19 vaccines approved by regulatory bodies around the world have been shown to provide significant protection against COVID-19. But no vaccine is 100% effective at preventing infection, also they have varying efficacy rates and different side effects. However, the main protease (Mpro) of SARS-CoV-2 has been identified as a key drug target due to its essential role in viral infection and its minimal similarity with human proteases. Cordyceps mushrooms have been found to have various therapeutic properties that could effectively combat SARS-CoV-2, including improve lung functioning, anti-viral, immunomodulators, anti-infectious, and anti-inflammatory. The present study aims to screen and evaluate the inhibitory potential of the bioactive molecules from the Cordyceps species against the Mpro of SARS-CoV-2. The bioactive molecules were screened based on their docking score, molecular interactions in the binding pocket, ADME properties, toxicity, carcinogenicity, and mutagenicity. Among all the molecules that were tested, cordycepic acid was the most effective and promising candidate, with a binding affinity of −8.10 kcal/mol against Mpro. The molecular dynamics (MD) simulation and free binding energy calculations revealed that the cordycepic acid-Mpro complex was highly stable and showed fewer conformational fluctuations. These findings need to be investigated further through in-vitro and in-vivo studies for additional validation. Communicated by Ramaswamy H. Sarma


Introduction
SARS-CoV-2 belongs to the family of large coronaviruses that commonly cause respiratory disease in humans and some animals (Zhang & Holmes, 2020).It bears a strong resemblance to the SARS-CoV that previously infected humans.The first case of SARS-CoV-2 was reported in Wuhan, China in late 2019 (Huang et al., 2020) and it rapidly spread to other major cities of China and other countries.In such an unprecedented situation with a high incidence of infections, the World Health Organization (WHO) declared the outbreak of this novel coronavirus a global pandemic on March 11, 2020.This declaration starts a new challenge for the development and production of new inhibitors against SARS-CoV-2 to end this pandemic.However, the development of an effective therapy, such as development of small molecules and the search for natural products that can be effectively used against SARS-CoV-2, is a challenging task for the scientific community.As of 5:43 PM on May 17, 2023, the worldwide cumulative count of confirmed COVID-19 cases has reached 766,440,796, accompanied by a total of 6,932,591 reported deaths.Consequently, the calculated mortality rate for COVID-19, based on these figures, is approximately 0.904%.
Genomic studies have shown that SARS-CoV-2 is closely related to other viruses, such as SARS-CoV and MERS-CoV, with a sequence similarity of 82.1% (Khan et al., 2020) which reaches up to > 90% for structural proteins and key enzymes (Ghosh et al., 2020).However, the drugs and other treatment used to fight against SARS-CoV and MERS-CoV were unsuccessful against SARS-CoV-2 (Naqvi et al., 2020).The viral genome is a positive single-stranded RNA that encodes for several structural and non-structural proteins (Low et al., 2021).The four major structural proteins are the spike protein (S), which is required for recognition and binding to the receptor, an envelope protein (E), which aids in assembly and maturation, membrane protein (M), which aids in assembly and budding, and nucleocapsid protein (N), which is required for replication and packaging.The spike protein of the mature virus consists of two non-covalently linked functional subunits, S1 and S2.The S1 subunit is responsible for binding to the angiotensin-converting enzyme-2 (ACE-2) receptor on host cells, which leads to the formation of endosomes and the fusion of viral particles.The S2 subunit, on the other hand, contains a fusion peptide that facilitates the anchoring of the virus to the host cell membrane, which is a critical step in the infection process, as shown in Figure 1.
Once a virus particle is internalized, its envelope is removed, and its genome is exposed to the host cell cytoplasm (Ruan et al., 2003).In the cytoplasm, positive-sense single-stranded RNA of the SARS-CoV-2 produced two large polyproteins, pp1a, and pp1ab, from the conserved ORF1a and ORF1ab genes as shown in Figure 2.These polyproteins are processed by two viral cysteine proteases: Papain-like protease (PLpro), which is responsible for proteolytic cleavage of nsp 1-4, and main protease (Mpro), which catalyzes proteolytic cleavage of nsp 5-16 (Zou et al., 2020).
SARS-CoV-2 Mpro, also known as 3CLpro or main protease, is an enzyme that plays a critical role in the replication of the virus.It is a homodimeric cysteine protease with a molecular weight of approximately 33.8 kDa per monomer.Each monomer of the Mpro consists of three domains: Domain I (residues 8-101), Domain II (residues 102-184), and  Domain .Domains I and II both include beta-barrels that have a chymotrypsin-like structure.These barrels also contain a catalytic dyad consisting of histidine 41 and cysteine 145 (Chan et al., 2020;Qiu & Xu, 2020).Domain III comprises five alpha helices and is involved in the dimerization of Mpro.This dimerization process is required for the synthesis of twelve nsp (Nsp5-Nsp16) (Fang et al., 2008;Kneller et al., 2020), including the most critical proteins, such as helicase and RNA-dependent RNA polymerase (RdRp).
Mpro play a very crucial role in the viral replication, as well as in the transcription of subgenomic RNA that required for the synthesis of viral structural proteins.Mpro is widely conserved among coronaviruses and exhibits a high degree of homology with other members of the same family such as MERS-CoV and SARS-CoV.It has been experimentally observed that inhibition of Mpro suppresses the infection rate of SARS-CoV-2 (Wong & Saier, 2021).For example, a bisulfate-based prodrug GC376 interact with Mpro and inhibits SARS-CoV-2 infections (Vuong et al., 2020).Similarly, boceprevir, a clinically approved drug of hepatitis C virus interact with Mpro and inhibits the infection of SARS-CoV-2 (Fu et al., 2020).It is also found that the virus may penetrate different host cells via distinct receptors.However, the activity of Mpro is required by the virus in all the cell types for proliferation and maturation (Perrier et al., 2019).Due to such significant involvement, Mpro has become a potential target for drug development against SARS-CoV-2.Several drugs that target the Mpro protease of SARS-CoV-2 like GC376, PF-07321332, AT-527, Ebselen, and Nafamostat are currently in clinical trials for the treatment of COVID-19.
Previous studies have shown that natural products could be a potential inhibitor of SARS-CoV-2.For example, the inhibitory effect of bioactive tea molecules against SARS-CoV-2 was investigated.In this study, it was reported that Oolonghomobisflavan-A is a potent inhibitor of the Mpro, based on molecular docking, molecular dynamics and free binding energy calculations (Bhardwaj et al., 2021).Similar study on bioactive tea molecules against Nsp15 SARS-CoV-2, reported that three molecules present in tea, namely barrigenol, kaempferol, and myricetin, are the potential lead molecules that showed the highest inhibitory potential against the non-structural protein-15 of SARS-CoV-2.(Sharma et al., 2021).Some previous reports also revealed that curcuminoids showed stable and high binding affinity towards Nsp15 SARS-CoV-2 and could be admitted as potential inhibitor (Singh et al., 2022).Various studies have highlighted the effectiveness of plant-derived molecules in inhibiting the infection of SARS-CoV-2 virus.These molecules have been found to exhibit strong antiviral properties that can potentially block the viral entry and replication mechanisms, thus preventing the onset of COVID-19.In addition to their potency, these plant-based inhibitors have also been recognized for their low toxicity levels, making them a safer alternative to synthetic drugs.Such findings indicate the potential for plant-based treatments to offer a viable solution in the fight against the ongoing COVID-19 pandemic (Khan et al., 2021(Khan et al., , 2022;;Trivedi et al., 2022).Among various natural sources, mushrooms have their own importance.From ancient time humans has used mushrooms for treating various diseases (Papaj et al., 2022).Cordyceps mushrooms, in particular, are reported to have immunomodulatory (Ng & Wang, 2005), anti-infectious (M€ uller et al., 1991), anti-inflammatory (Kim et al., 2006), anti-tumor (Nakamura et al., 2006), andanti-viral properties (WenJuan et al., 2010).Cordyceps has significant role in the treatment of respiratory diseases like chronic obstructive pulmonary disease, bronchitis and asthma.It has also been reported to induce proliferation of T cells, cytotoxic activity of natural killer (NK) cells, activates macrophages and secretion of Th1 cytokines.(Wook Baik Bundang Jesaeng Hospital et al., 2015).One of the previous studies reported that how cordycepin (3 0 -deoxyadenosine), a bioactive molecule of the Cordyceps species can be used as an activator for the adenosine receptors that helps to enhance immunity against SARS-CoV-2 (Du et al., 2021).
In the era of this pandemic, Cordyceps may be a key therapeutic option to reduce the infection of SARS-CoV-2.The present study aims to screen and evaluate the inhibitory potential of the bioactive molecules from the Cordyceps species against the Mpro of SARS-CoV-2.The first objective of the study is to analyze the binding affinity of the molecules with Mpro using molecular docking.Second, to evaluate drug-likeness properties based on the Lipinski rule of five of the screened molecules.Third, to analyze the toxicological, carcinogenic, and mutagenic properties of the studied molecules.Finally, the in-depth analysis of protein-ligand complexes using molecular dynamics simulations and binding free energy calculations.

Protein structure and its preparation
The X-ray crystal structure having resolution 1.59 Å of the Mpro (PDB ID 5R7Z) (Douangamath et al., 2020) was retrieved from the PDB (Protein Data Bank).The protein preparation wizard of the Schr€ odinger suite was used in default mode to incorporate zero bond order for metal, assign proper bond orders, generate disulphide bonds and add missing hydrogen atoms to the target protein.All water molecules are removed in the radius of 3 Å from HET groups.Missing loops and side chains were added using the prime module of the Schr€ odinger suite.PROPKA program used for evaluate the pKa of protein residues at pH 7.0.The target protein was further optimized to form an H-bond network.The OPLS3e force field was employed to minimize the energy with the RMSD of 0.30 Å.

Ligand preparation
The 3-dimensional structures of the bioactive molecules were retrieved from the PubChem database.All the ligands were prepared using LigPrep (Schr€ odinger Release 2022-2).Epik program of Schr€ odinger suit was used to generate ionization states of the ligands at pH 7.0 ± 2.0.A maximum of 32 stereoisomers were generated for each ligand, and energy minimization was carried out using the OPLS3e force field (Table 1).

Grid generation
The 'Glide's Receptor Grid Generation' program of the Schrodinger suit was used to create the grid box around the binding site of the co-crystallized native ligand with the dimensions of X: 11.38 Y: 2.69 Z: 23.14 with van der Waals scaling factor of 1.0 and partial charge cutoff at 0.25.

Molecular docking
The ligand set, which included co-crystallized ligand, was prepared and subjected to molecular docking against Mpro using the extra precision mode (XP) with flexible ligand sampling in the Maestro Glide docking module (Friesner et al., 2006).The docking utilized the OPLS3e force field, with a van der Waals scaling factor of 0.80 and a partial charge cutoff at 0.15.The molecules were assessed based on their docking scores, indicating their interaction affinities, and subsequently chosen for additional analysis.

Assessment of ADME properties
e SwissADMEdatabase and DruLiTo software was utilized (Daina et al., 2017) to calculate the ADME properties of the compounds, along with their associated properties.To evaluate the drug-likeness of the compounds, the Lipinski rule of five was applied, which is based on molecular weight, LogP value, hydrogen bond donors, and hydrogen bond acceptors.

Assessment toxicity, carcinogenicity and mutagenicity using QSAR
In this study, pkCSM database was employed to evaluate the toxicity of the screened molecules (based on molecular docking and ADME assessment) (Pires et al., 2015).To predict the carcinogenicity and mutagenicity of the compounds, Toxtree 3.1.0software was used to encode the Benigni/Bossa rulebase (Benigni et al., 2008).

Molecular dynamic (MD) studies
In the present study, Chemistry at Harvard Macromolecular Mechanics 27 (CHARMM27) was used for molecular simulation program to investigate the molecular dynamics of the docked complex.CHARMM27 is a widely used and highly versatile tool (Brooks et al., 2009).Additionally, e Nanoscale Molecular Dynamics (NAMD) version 2.8 (Nelson et al., 1996) and the Transferable intermolecular potential with 3 points (TIP3P) model for water (Jorgensen et al., 1983) was employed for-.The topology and parameters of the ligand are generated using SwissParam (Zoete et al., 2011), and TIP3P water molecules were added to fully cover the entire complex surface, which had approximate dimensions of 98.5 � 98.5 � 69.7.To neutralize the net charges of the system, we added sodium (4) and chloride (0) counterions.
After setting up the solvent environment, each complex system consisted of approximately 67,300 atoms.Molecular dynamics (MD) simulations for 100 ns, first equilibrating the equations of motion for 50,000 steps (100 ps) in a constant temperature and constant volume (NVT) ensemble at 300 K, followed by the NPT ensemble for an additional 50,000 steps (100 ps).During the production dynamics phase, the isobaricisothermal, constant temperature and constant pressure (NPT) ensemble with a 2fs time step and saved coordinates every picosecondwere used for analysis.To visualize the protein-ligand complex and analyze the MD trajectory, Visual Molecular Dynamics (VMD) version 1.9.1 software was employed (Humphrey et al., 1996).

Binding free energy (Dg bind ) studies
The Molecular Mechanics-Generalized Born Surface Area (MMGBSA) and Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) methods, both incorporated in the MMPBSA.pymodule of AMBER18 (Miller et al., 2012), were utilized to estimate the free energy of binding for the docked complex.After analyzing 100 frames from the trajectories, the net energy of the system was calculated using the following equation:

Molecular docking
Study utilized the GLIDE module of Schr€ odinger to perform docking studies with a prepared target protein, 5R7Z, and a dataset comprising bioactive molecules from Cordyceps species.In order to set up the docking grid, we centered it around the co-crystallized ligand, which had the following coordinates: X: 11.38, Y: 2.69, Z: 23.14.To assess the potential binding affinity of each molecule, a docking score was assigned, and only those with a score of À 5.0 kcal/mol or lower were considered appropriate for further analysis (Table 2).
Ensuring the reliability and accuracy of molecular docking simulations is crucial in drug design, and this necessitates the validation of the docking protocol.Validation involves comparing predicted binding modes and affinities with experimental data or reference structures to identify potential sources of error and optimize the docking protocol parameters to enhance its performance.Furthermore, validation can help to identify false-positive or false-negative results and determine the sensitivity and specificity of the protocol.In the present study, the docking protocol was validated by extracting the native cocrystalized ligand from 5R7Z and redocked at the binding site, as illustrated in Figure 3.The position, conformation, and orientation of the redocked cocrystlized ligand were found to be similar to the native structure obtained from the PDB.The accuracy of the docking protocol was confirmed by the generation of the same bonding interactions and an RMSD value of 0.7.

Interaction studies
Among the data sets that were docked, cordycepic acid displayed the best binding position, with a docking score of À 8.100 kcal/mol, followed by cicadapeptin II and N-2-hydroxyethyl adenosine with docking scores of À 7.149 kcal/mol and À 6.941 kcal/mol, respectively.Cordycepic acid binds tightly to five hydrogen bonds with residues GLU166 (2H-B), ARG188, GLN189, and THR190.Additionally, it interacts with residues MET165, GLU166, LEU167, PRO168, GLN192, ALA 191, THR190, GLN189, and ARG188 through hydrophobic, hydrophilic, and van der Waals interactions.Cicadapeptin II interacts with the binding pocket of Mpro through five hydrogen bonds with residues ASN142, GLU166, THR190, and GLN189.Furthermore, it forms a salt bridge with residue GLU166.N-2-hydroxyethyl adenosine interacts with the active site of Mpro through four hydrogen bonds with residues HIE41, GLY143, GLU166, and ARG188.N-2-hydroxyethyl adenosine also interacts with residue HIE41 through Pi-Pi stacking.A Cys145-His41 catalytic dyad facilitates the processing of viral polyproteins by Mpro.Of all the molecules, ten interacted with the residues of the catalytic dyad and other critical amino acids in the catalytic site.The binding site residues involved in the interactions between Mpro and bioactive molecules are comprehensively detailed in Figure 4A-J and Table 3.

Assessment of ADME properties
ADME properties play a crucial role in drug development, as they determine how a drug is absorbed, distributed, metabolized, and excreted in the body.These properties can affect a drug's efficacy and safety in humans.For example, if a drug is not absorbed efficiently, it may not reach its target site in the body, rendering it ineffective.Similarly, if a drug is metabolized too quickly or excreted too rapidly, it may not reach a sufficient concentration in the body to produce the desired effect, or it may accumulate to toxic levels.Therefore, evaluating ADME properties is critical in determining the optimal dose, route of administration, and dosing frequency of a drug, as well as identifying potential safety concerns before its approval for human use.This information can help to ensure that the drug is safe, effective, and has a higher chance of success in clinical trials.Overall, the assessment of ADME properties is a crucial step in drug development that helps ensure the safety and efficacy of new drugs.The Lipinski Rule of 5 is a tool used to predict the likelihood of a compound becoming a successful drug.
Compounds that violate the rule are less likely to have good pharmacokinetic properties and may have poor absorption, distribution, metabolism, and excretion (ADME) properties.This rule includes the number of hydrogen bond acceptors (� 10), the number of hydrogen bond donors (� 5), molecular weight (� 500), and the log Po/w (� 5).A molecule with more than one rule violation is not allowed in molecular dynamic studies.Therefore, Lipinski Rule of 5 ADME analysis is an essential part of drug discovery, as it helps identify compounds that are more likely to be developed into successful drugs, potentially saving time and resources.Lipinski's rule of five can help to determine whether a biologically active molecule will have the chemical and physical properties necessary for oral route.Among ten molecules analyzed, two molecules, Cicadapeptin II and Ophiocordin violated more than three rules.Log Po/w is a crucial parameter that help to understand the pharmacokinetic properties of the molecule and dictates how a drug should be formulated and dosed.Six molecules, Adenosine, Cordycepic acid, Cordycepin, Guanosine, N-2-Hydroxyethyl-adenosine and Uridine have negative log Po/w values, which indicates they have a higher affinity for the aqueous phase.Such molecules are readily soluble in the blood, which increases their bioavailability.On the other hand, molecules with low aqueous solubility (high positive log Po/w) will be compromised in bioavailability; also, they will not be able to cross the blood-brain barrier.Assessment of ADME properties of all the molecules are summarized in the Table 4.

Assessment of toxicity, carcinogenicity and mutagenicity using QSAR
Assessing the potential toxicity, carcinogenicity, and mutagenicity of new drugs is crucial to ensure their safety and efficacy.Toxicity refers to the harmful effects a drug can have on living organisms, including humans.Carcinogenicity refers to the drug's potential to cause cancer, while mutagenicity refers to its potential to cause genetic mutations.Identifying these potential effects during the drug development process is essential to ensure the drug's safety for human use.This information can help in making informed decisions about the drug's development, including whether to continue development, modify its structure or dosages, or halt its development altogether.Overall, assessing toxicity, carcinogenicity, and mutagenicity is critical for ensuring the safety and effectiveness of new drugs and protecting public health.To evaluate the toxic, carcinogenic, and mutagenic potential of the screened molecules, the pkCSM database in toxicity mode and Toxtree 3.1.0,using the Benigni-Bossa rulebase, were employed.Among the ten molecules studied, N-2-hydroxyethyl-adenosine, cicadapeptin II, cordycepin, and guanosine were found to be hepatotoxic, while myriocin, adenosine, cicadapeptin II, cordycepin, and guanosine were identified as genotoxic carcinogens.Additionally, cicadapeptin II was classified as a non-genotoxic carcinogen.Four molecules from the cordyceps species, cordycepic acid, lovastatin, ophiocordin, and uridine, were analyzed and found to be non-toxic, non-carcinogenic, and non-mutagenic, as illustrated in Figure 5 and Table 5.

Molecular dynamics and binding free energy (Dg bind )
Cordycepic acid was identified as the best lead molecule among all screened bioactive molecules from Cordyceps species based on docking scores, interaction studies, ADMET properties, and carcinogenic-mutagenic profiling.To further evaluate the dynamic behavior and stability of the Cordycepic acid and Mpro docked complex, a 100 ns Molecular dynamics (MD) simulation was conducted.The MD trajectory was used to calculate various parameters such as root mean square deviation (RMSD), root mean square Cicadapeptin II À 7.149 3.
The RMSD value is a significant parameter used in molecular dynamics simulations to evaluate the stability of a protein-ligand complex.It determines the distance between the atoms of a simulated structure and a reference structure, like a crystal structure.A low RMSD value suggests a stable interaction between the molecules since the simulated structure is in close agreement with the reference structure.Conversely, a high RMSD value suggests an unstable interaction due to structural deviation from the reference structure.The RMSD value is valuable in assessing the quality of a molecular dynamics simulation and optimizing the simulation parameters to obtain a stable structure close to the reference structure.In this study, the stability of the protein-ligand interaction was validated using the RMSD value.The simulation showed an average backbone RMSD of 0.21 nm, ranging from 0.04 to 0.29 nm.However, fluctuations were observed throughout the production run, as depicted in Figure 6.
The root-mean-square fluctuation (RMSF) is an important parameter used in molecular dynamics simulations to measure the mobility and flexibility of molecules.It measures the average displacement of each atom in a molecule relative to its average position during simulation and can be used to identify flexible or rigid regions of a molecule.RMSF values can also be used to monitor the structural stability of proteins and other biomolecules, and compare flexibility between different molecules.In this study the complex slightly fluctuated between the residue 2400-2500 and 4200-4300.The fluctuation in the terminal part is higher than that in the central part, as shown in Figure 7.Such behavior was expected because the central region, which has a more rigid alpha-helix and beta-strand structures, fluctuates less than the terminal region, where fluctuation is always at a higher level due to the free ends.In the course of the simulation, an average RMSF of 0.09 nm was observed, with a range spanning from 0.03 to 0.55 nm.According to the RMSF findings, the complex is stable with low conformational flexibility.
The radius of gyration is a critical parameter in molecular dynamics simulations that provides valuable information about the compactness and conformational changes of a molecule.It is calculated as the root-mean-square distance of all atoms in a molecule from their center of mass, and it reflects the size and shape of the molecule.The radius of gyration can be used to monitor protein folding and unfolding during simulation, as a decrease indicates a more compact and folded conformation, while an increase indicates an unfolded or extended conformation.In this simulation study of 100 ns, the radius of gyration values started at �2.14 nm and increased up to �2.19 nm, indicating the stability and integrity of the complex.However, the proline at the 168th position acted as a hinge or separator of the two beta sheets (Figure 8).
Hydrogen bonds are essential for maintaining the stability and structure of biomolecules during molecular dynamics simulations.They can influence conformational changes, which can alter the biomolecule's function.Monitoring hydrogen bond formation and breaking provides valuable insights into the dynamics of a system.If hydrogen bonds are stable, the protein or biomolecule maintains its structure, while the breaking of hydrogen bonds may indicate unfolding or conformational changes.In protein-ligand interactions, hydrogen bonds are critical for binding affinity and specificity.The formation of hydrogen bonds between a ligand and a protein's active site enhances binding, while their breaking can reduce binding affinity.The study employed a cutoff value for distance and angle of 3.5 Å and 30 � , respectively, to calculate the number of hydrogen bonds formed between a ligand and protein during a 100 ns MD simulation.The results showed that the ligand had a strong affinity for the target protein, with a maximum of 6 hydrogen bonds observed and an average of 2 hydrogen bonds formed between the ligand and protein.The findings suggest that the ligand is a potential candidate for drug development against the target protein (Figure 9).
The Solvent Accessible Surface Area (SASA) method is a popular technique utilized in molecular dynamics simulations to assess conformational changes in a protein and its interactions with solvent molecules.It computes the protein's surface area that can be accessed by the solvent molecules, providing valuable insights into the protein's dynamics and stability in solution.In the current study, SASA was employed to evaluate the conformational changes in the protein structure for 100 ns in MD simulation, as demonstrated in Figure 10.The average SASA value was calculated to be 149.2nm 2 , indicating that the complex remained stable throughout the simulation.This suggests that the SASA method is a reliable tool for analyzing protein behavior in solution and can have significant implications for protein-ligand interactions and drug design.
MMGBSA (Molecular Mechanics-Generalized Born Surface Area) and MMPBSA (Molecular Mechanics-Poisson Boltzmann Surface Area) are widely used methods to estimate the binding free energy of ligand-protein complexes.
These methods consider various energy components, including van der Waals, electrostatic, and solvation energies, to provide a comprehensive understanding of the binding interactions.The binding energy, which refers to the energy released during the formation of chemical bonds or the interaction between a ligand (a molecule that binds to a specific site on a protein) and a protein.The strength of this interaction is measured by the binding energy, with a lower binding energy indicating a better interaction between the ligand and protein.In order to verify the ligand affinity for the receptor protein complex predicted by the docking simulation studies, the binding free energy of the complex was calculated.The findings of the study revealed that electrostatic interactions dominated the binding energy, with a value of À 44.616 kcal/mol, followed by van der Waals interactions with a value of À 26.254 kcal/mol.The polar solvation energy, which represents the energy required to dissolve polar molecules in a solvent, was unfavorable in the formation of the cordycepic acid and Mpro complex, with values of 26.580 kcal/mol and 31.276kcal/mol for MMGBSA and MMPBSA, respectively.However, despite the unfavorable polar solvation energy values, the complex predicted by MMGBSA had a net binding energy of À 30.070 kcal/mol, indicating greater system stability than the complex predicted by MMPBSA, which had a net energy of À 26.708 kcal/mol.The results of the study, including the binding free energies of the system for MMGBSA and MMPBSA, are summarized in Table 6.

Conclusion
This study screened and evaluated bioactive molecules from Cordyceps species for their inhibitory potential against the SARS-CoV-2 main protease.The binding pattern of these molecules with Mpro was investigated using molecular docking.Several bioactive molecules, including cordycepic acid, cicadapeptin II, N-2-hydroxyethyl adenosine, myriocin, uridine, ophiocordin, guanosine, adenosine, cordycepin, and lovastatin, were found to have strong binding with Mpro.These molecules were further evaluated based on ADMET properties and carcinogenic-mutagenic profiling.
Cordycepic acid was found to be the most promising candidate and was subjected to molecular dynamics simulations.The results of the MD simulations, including RMSD, RMSF, Rg, number of H-bonds, and SASA, suggested that the docked complex was highly stable.Furthermore, MMGBSA and MMPBSA calculations of binding free energy were used to validate the results and evaluate the binding affinity of cordycepic acid to Mpro.
The study concludes that cordycepic acid is a novel inhibitor of SARS-CoV-2 main protease.However, further in-vivo      and in-vitro analyses are required to investigate its inhibitory potential against Mpro.Overall, the present study highlights the potential of natural products, such as those derived from Cordyceps species, as a source of potential inhibitors against SARS-CoV-2 main protease.

Table 3 .
Interactions of bioactive molecules and Mpro.

Table 4 .
Assessment of ADME properties.Molecule with more than one violation of the rule is not allowed in the molecular dynamic studies.

Table 6 .
Binding free energies for Cordycepic acid and Mpro complex.