Computational approach investigation bioactive molecules from Saussurea Costus plant as SARS-CoV-2 main protease inhibitors using reverse docking, molecular dynamics simulation, and pharmacokinetic ADMET parameters

SARS-COV-2 virus causes (COVID-19) disease; it has become a global pandemic since 2019 and has negatively affected all aspects of human life. Scientists have made great efforts to find a reliable cure, vaccine, or treatment for this emerging disease. Efforts have been directed towards using medicinal plants as alternative medicines, as the active chemical compounds in them have been discovered as potential antiviral or anti-inflammatory agents. In this research, the potential of Saussurea costus (S. Costus) or QUST Al Hindi chemical consistent as potential antiviral agents was investigated by using computational methods such as Reverse Docking, ADMET, and Molecular Dynamics with different proteases COVID-19 such as PDB: 2GZ9; 6LU7; 7AOL, 6Y2E, 6Y84. The results of Reverse Docking the complex between 6LU7 proteases and Cynaropicrin compound being the best complex, as the same result, is achieved by molecular dynamics. Also, the toxicity testing result from ADMET method proved that the complex is the least toxic and the safest possible drug. In addition, 6LU7-Cynaropicrin complex obeyed Lipinski rule; it formed ≤5 H-bond donors and ≤10 H bond acceptors, MW < 500 Daltons, and octanol/water partition coefficient <5.


Introduction
Since (2019) COVID-19 virus has been announced as a global pandemic, COVID-19 by the World Health Organization (W.H.O), on February 11th, 2020 stands for "Co": corona (meaning crown), "vi": virus, "d": disease, and "19": the year (2019) [1][2][3][4][5]. This virus which is labeled as (COV2) is different from both MERS-COV (Middle East Respiratory Syndrome-related Coronavirus) responsible for an epidemic that spread in 2012 in the Middle East [6], and SARS (Severe Acute Respiratory Syndrome Corona Virus1 that brought SARS epidemic about in 2003 [7]. The virus enters the body through contact with the cells of the mucous membranes lining the nose, mouth, and eyes, causing infections in the respiratory system, and it takes at least six days (the incubation period) [8,9] for symptoms to appear in the infected patients.
The main symptoms of COVID-19 disease are extreme fatigue, fever, sneezing, shortness of breath, coughing, and breathing difficulties. These symptoms are similar to those of a cold or the flu. Still, the impact of COVID-19 is different from one person to another, as the risks of severe complications of COVID-19 that led to death or prolonged recovery with weak lungs in smokers and elderly patients [10]. Zhang et al. examined the mechanism of SARS-CoV-2 virus action with human body, it exploits the interplay between the small miRNA and other biomolecules to avoid being effectively recognized and attacked by host immune protection and deactivate functional genes crucial for the immune system. In detail, SARS-CoV-2 can be regarded as a sponge to adsorb host immune-related miRNA, which forces the host to fall into the dysfunctional status of the immune system. Besides, SARS-CoV-2 encodes its miRNAs, which can enter a host cell and are not perceived by the host's immune system, subsequently targeting host function genes to cause illnesses. Therefore, this article presents a reasonable view that the miRNA-based interplays between the host and SARS-CoV-2 may be the primary cause that SARS-CoV-2 accesses and attacks the host cells [35,36].
During the epidemic, there was an increase in the consumption of Saussurea costus in Indonesia, which has a low rate of hospital admissions for Indonesian patients compared to other countries [11]. The Costus plant originated in India and has since spread to many countries worldwide. It has a ginger-like appearance, grows to a height of 1 m and a half, and has been used for medicinal purposes since antiquity [12].
There are two types of Costus: the white, sweet "Marin" or "Bahri," which the Arabs brought to their countries by sea and whose taste is prized, and the black, bitter "Indian," which is thought to be superior in terms of medicine strength and defense against a wide range of diseases, especially influenza [13]. The most important part used in the treatment is the skin of the roots [14,15].
Costus powder is usually used for medical purposes, where the active ingredients are organic compounds. It has a pivotal role in treating allergies and is an antiseptic agent for germs. It is an essential antiseptic   agent for tonsillitis germs, pharyngitis, bacterial pleurisy, colds, influenza, pneumonia, and several germs that infect the lungs. In addition, it has a significant role in shortness of breath, asthma, and cough. Because it prevents and mitigates harmful and sometimes fatal attacks of viruses on the respiratory system and supports the immune system. The oil and roots of Costus are the other parts of the plant used for medical purposes, which must be used scientifically and pharmacologically to avoid overdose poisoning [16].
This research aims to investigate Costus plant active ingredients as a potential treatment for Covid-19 relying on the pharmaceutical properties of Saussurea Costus. The current study uses different computational approaches; such as Reverse docking, Molecular dynamics, and ADMET. In addition, Lipinski rules to screen the interest of this medicinal plant against SARS-CoV2 infection (see Fig. 1).

The bioactive molecules of the medicinal plant Saussurea costus
Three organic bioactive compounds have led to increased interest in researchers, as shown in Fig. 2 such as Costunolide, Dehydrocostus lactone, and Cynaropicrin are sesquiterpene lactones found in Saussurea costus contains the most important amino acids, organic acids, essential oil, highlighting strong potential therapeutic effects helps to recover faster with various acute respiratory diseases, flu, sore throat and pleurisy [18][19][20].

Reverse docking
Recently, the development of computer and bioinformatics programs has made great strides in predicting the molecular interactions that carry a protein and ligand at the binding site and in designing and discovering new drugs. The bioinformatics program that we applied in this research is reverse docking molecular modeling. This program is widely used mainly in medicinal chemistry and in designing new therapeutic molecules [21]. By fixing the target molecules at the disease protein's binding sites and looking for the most stable interactions, reverse docking aids in the discovery of the best position for the bioactive molecule and disease protein. The molecular structures of Costunolide, Dehydrocostus and Cynaropicrin were energetically minimized and geometrically optimized by using Gaussian software. Then the realization of reverse docking is also based on the preparation of proteins: The crystal structures of SARS-CoV-2 are available in the Protein Data Bank (PDB, htt ps://www.rcsb.org/), the five models for high-resolution SARS-CoV-2 (PDB ID: 2GZ9 (2.17 Å); 6LU7 (2.16 Å); 7AOL (1.47 Å), 6Y2E (1.75 Å), 6Y84 (1 0.39 Å) were downloaded from the Protein Data Bank Fig. 3. These proteins were prepared by removing water molecules and adding polar hydrogens and Kohlman charges by using the Discovery Studio program, the prepared files are then converted to PDB format. As well as the sybyl surflex-dock software is used for the calculations and the study of molecular docking. The studied complexes were analyzed and reassembled in PDB format by using PyMol software [22,23]. Finally, they were sent to the Discovery Studio software to visualize the interactions [24].

Molecular dynamics simulations
GROMACS simulation package (GROMACS 2020.4) used to perform molecular dynamics simulations for protein-ligand complex was carried out for 100 ns in water using CHARMM36 force field; trajectory and energy files written every 10 ps [25].
The system was solvated in a truncated octahedral box containing TIP3P water molecules. The protein was centered in the simulation box within a minimum distance to the box edge of 1 nm to satisfy the minimum image convention efficiently. The simulation was performed in 0.15 M KCl by adding 60 Potassium ions and 57 Chloride ions, and the overall system contained 66,242 atoms.
Minimization was carried out for 5000 steps using Steepest Descent Method, and the convergence was achieved within the maximum force <1000 (KJ mol − 1 nm − 1 ) to remove any steric clashes. All three systems were equilibrated at NVT and NPT ensembles for 100ps (50,000 steps) and 1000ps (1,000,000 steps), respectively, using time steps 0.2 and 0.1 fs, at 300 K to ensure a fully converged system for a production run.
The production runs for simulation were carried out at a constant temperature of 300 K and a pressure of 1 atm or bar (NPT) using weak coupling velocity-rescaling (modified Berendsen thermostat) Parrinello-Rahman algorithms, respectively. Relaxation times were set to τ T = 0.1 ps and τ P = 2.0 ps. All bond lengths involving hydrogen atoms were kept rigid at ideal bond lengths using the Linear Constraint Solver (lincs) algorithm, allowing for a time step of 2 fs. Verlet scheme was used for the calculation of non-bonded interactions. Periodic Boundary Conditions (PBC) were used in all x, y, and z directions. Interactions within a short-range cutoff of 1.2 nm were calculated in each time step. Particle Mesh Ewald (PME) was used to calculate the electrostatic interactions and forces to account for a homogeneous medium outside the long-range cutoff. The production was run for 100 ns for the complex [26].

Predicted pharmacokinetic and toxicity properties
It is the science of drug kinetics Absorption, Distribution, metabolism, and elimination (kinetic study of ADMET), [A]bsorption the passage of the drug from the site of administration into the general circulation.
[D]istribution is the movement of blood to and from the tissues.
[M]etabolism of drugs is the chemical modification of the drug with the overall goal of getting rid of the drug; enzymes are generally involved in metabolism.
[E]limination is the irreversible removal of the parent drug from the body.
[T]oxicity as the name suggests, this filter is used to measure the value of a compound and its metabolites [27][28][29].
Lipinski's rules are a list of guidelines proposed by Chris Lipinski of Pfizer to help design molecules with good availability. The four guidelines are: (1) The molecule must have a molecular weight MW of 500 g/ Mol or less.
(2) the molecule should have no more than 5 hydrogen bond donors HBDs. These are mainly OH and NH groups. (3) Molecules should have no more than 10 HBAs than a hydrogen bond acceptor. They are essentially the oxygen and nitrogen atoms in the molecule (4). The lipophilicity of the molecule, measured by logP, must not be greater than 5. These have been called "Lipinski's Rules" or "Rule of 5" [30,31]. Lipinski recommended that the molecules adhere to at least three of the four guidelines to have a reasonable chance of having beneficial oral bioavailability. Therefore, when we think about ADME, Lipinski's rules help meet the letter A-absorption or have higher efficacy or lower toxicity. Therefore, these rules focus only on improving solubility and membrane permeability to maximize oral absorption.  Computational approaches can help minimize the risks of toxicity. In this research, a new approach pkCSM has developed. It is considerably better than current methods, where the freely accessible web servers (SwissADME software) that Simplified Molecular Input Line System (SMILES) used during the generation process [32].

Reverse docking
Reverse docking was applied on three molecules of Saussurea costus Costunolide, Cynaropicrin, and Dehydrocostus Lactone to identify therapeutic activity against Sarc-Cov-2 [18]. It is a technique that allows predicting or studying the interactions that are the main factors that significantly impact the affinity of a ligand for a receptor. In addition, this study helps us choose the effective protein for some compounds based on Saussurea costus, considered a treatment against COVID-19. The five proteins that are available in the Protein Data Bank with good resolution (ID PDB: 2GZ9 (2,17 Å); 6LU7 (2,16 Å); 7AOL (1,47 Å), 6Y2E (1,75 Å), 6Y84 (1,39 Å) on three molecules of Saussurea costus Costunolide, Cynaropicrin, and Dehydrocostus as shown in Fig. 4. For PDB 2GZ9 (2.17 Å) (Fig. 5), it was observed that there are just Pi-Alkyl-type interactions with the residues PHE 8 and PHE 294 in the case of Costunolide. On the other hand, visualization of docking of Cynaropicrin docked with the same protein gives best binding modes of protein with green color for Amino acid residues involved in hydrogen bonding interaction with residues THR 111 and TYR 154, As well, this complex provides another interaction Pi-Alkyl with residue PHE 294. Moreover, Dehydrocostus lactone shows Carbon-Hydrogen Bond interactions with residues PRO 108, Pi-Alkyl with ILE 200, LEU 202, and HIS 246. In addition, the PDB: 6LU7 complex (Resolution: 2.16 Å) (Fig. 6) with Cynaropicrin causes hydrogen bonds with the OH groups, which shows importance of these amino acids THR 111, ASP 153, and  SER 158, which are considered among the essential factors of stability protein-ligand. Moreover, the distances of its bonds are lower than 4 Å, which is confirmed by the strong hydrogen bonds; they are 1.56 Å, 2 Å and 3.34 Å respectively. Regarding 6LU7 protein with these two molecules, Costunolide and Dehydrocostus Lactone made Pi-Alkyl interactions with phenylalanine. On the one hand, the presence of hydrogen bond interaction of complex 7AOL (Resolution: 1.47 Å) (Fig. 7) and Costunolide with residue LYS 97 at a distance of 2.60 Å, a Carbon Hydrogen Bond interaction with ALA 70 at a distance of 2.52 Å and two Alkyl interactions with TRP 31 and VAL18 at distances of 4.46 Å and 5.01 Å, respectively. Thus, with Cynaropicrin, there are two interactions: Conventional Hydrogen Bond and Carbon-Hydrogen Bond with the amino acids GLN 110, SER and 158 successively.
Dehydrocostus Lactone is associated with the presence of two interactions formed with the amino acids of 7AOL, Pi-Alkyl with PHE 8 at a distance of 4.86 Å, and Pi-Donor Hydrogen Bond with PHE 294 at distances 3.26 Å.
On the other hand, the reverse docking study for the rest of the proteins, such as 6Y2E (Resolution 1.75 Å) and 6Y84 (Resolution 1.39 Å) (Figs. 8 and 9), give unsatisfactory results for Cynaropicrin. Also, it showed the presence of unfavorable donor-donor interactions with the amino acids LYS 102 and PHE 294, respectively.
According to molecular docking studies, the highest total score strongly denotes strong binding affinity. On the other hand, the molecule which will have minimum energy is considered the good one, so according to the two Tables 1 and 2, the best complex is [6LU7-Cynaropicrin] because it has a high score of 4.78 and low energy compared to the others with a value of − 7.4 kcal/mol and is not the only criteria, as well the best complex is [6LU7-Cynaropicrin] because of its chemical structure and docking interactions, as shown in (Fig. 10) is a guaianolide type known as a sesquiterpene lactone. It has a 5-7-5 fused tricyclic skeleton with four exo-olefins, six stereocenters, and two hydroxyl groups, which are caused by meaningful interactions of conventional Hydrogen Bond with THR 111, ASP 153 residues. Also, Pi-Alkyl and Alkyl interact with the amino acids PHE294 and ILE 106. Using molecular dynamics simulation, the structure of the SARS-CoV-2 main protease (PDB ID 6LU7) complex with the inhibitor cynaropicrin to confirm its stability.

Validation of docking protocol
The molecular docking validation step is crucial. In this research, the re-docking method was used, where a reference ligand is downloaded with the target protein to PYMOL software, which forms a new complex (docked) [37][38][39][40]. Fig. 11 demonstrates that the two ligands are positioned identically. In this operation, the small RMSD value of 1.986 Å indicates the accuracy of the docking results, which is close to the crystallography results.

RMSD (protein)
RMSD calculated the complex based on 'C-alpha' atoms using gromacs program. As shown in (Fig. 12), the C-terminal of the protein appears to be highly flexible (as indicated in RMSF), due to which the RMSD shows some fluctuation within a range of ~1.5 Å (in black). RMSD of protein (in Red) after removing the last 6 residues of C-terminal to check their effects on RMSD graph. The RMSD plot in red shows a relatively stable protein. The mean RMSD value for protein in black is 2.19 ± 0.43 Å and for protein in Red is 1.91 ± 0.15 Å.

RMSD (ligand)
RMSD is calculated for the ligand-based on ligand's atoms using gromacs program. The mean RMSD value is 1.25 ± 0.42 Å (see Fig. 13). RMSD of the ligand shows high fluctuations of the ligand. It also indicates that ligand exists most of the time in free form and not in bound form.

RMSF
A RMSF was calculated for the protein complex based on 'C-alpha' atoms using gromacs program. The whole protein remained relatively stable except for the C-terminal, which showed significant fluctuation (see Fig. 14).

ROG
Radius of gyration was calculated for 6LU7-Cynaropicrin (proteinligand) complex based on 'C-alpha' atoms using the gromacs program. The mean RoG value is 22.58 ± 0.22 Å (see Fig. 15). A very slight fluctuation is observed, indicating a slight opening and closing of protein. C-terminal plays a significant role in this fluctuation. Overall, the structure remains stable.

Protein-ligand complex
100 snapshots of ligand (Cynaropicrin) taken at every 1ns during 100 ns of simulation time. The image indicates that the ligand does not bind with the protein (6LU7). The initial conformation of (Cynaropicrin) ligand (0ns) is shown as red licorice. The conformation of (6LU7) protein is set to the initial conformation referring to 0ns of simulation (see Fig. 16). Fig. 17 shows the number of hydrogen bonds between 6LU7 and Cynaropicrin (a protein and linker) during 100 ns of simulation time. Both 6LU7 and Cynaropicrin show a slightly decreased number of hydrogen bonds during the simulation. However, Cynaropicrin is still in the whole mimic-related mode, which presents a link to other interactions that stabilize Cynaropicrin by the 6LU7 protein. Fig. 18 indicates that the Cynaropicrin (ligand) leaves the binding site and roams around the protein, with multiple short interactions. The more considerable COM distance continues until 90 ns, and later it comes closer to 6LU7 (protein).

Principal Component analysis (PCA)
The PCA of the formed 6LU7-Cynaropicrin complex calculated from the Bio3D application of R. is shown in Fig. 19. All three PCs captured 56.6% of the structural variance in protein. Fig. 20 shows atoms are  colored on a scale from blue to red, where blue represents atoms showing large motion amplitudes, and red represents more rigid atoms. The(c) shows the overlay conformations of 6LU7 (proteins) along both PC1 (shown as red) and PC2 (shown as blue), which are orthogonal to each other.

Dynamic cross correlation matrix analysis (DCCM)
As shown in Fig. 21 the colors varying from red to white to blue indicate the intensity of correlated motion. Blue colors indicate a negative correlation, white shows no correlation, and pink shows positively correlated motions between residues.

Potential energy
The 6LU7-Cynaropicrin complex system potential energy during 100 ns of MD simulation as obtained from Gromacs energy file. Fig. 22 shows the converged potential energy.

Pressure
The 6LU7-Cynaropicrin complex pressure and temperature at a 100 ns point in the MD simulation, as recovered from the GROMACS file, are presented in (Figs. 23 and 24). At some point throughout the 100 ns simulations, the diagram shows converged pressure and temperature.

MMGBSA binding energy
The binding energy of ligand-protein complexes during 100 ns of MD simulation was obtained by the MMGBSA method, where the value was − 1.74 ± 3.92 kcal/mol (see Fig. 25 & Table 3). The binding free energy is a crucial factor in determining the activity of drug molecules, and the lower value of ΔG bind value, the more stable 6LU7-Cynaropicrin complex formed. The MM/GBSA computation was once implemented by the use of MolAICal software. The Molecular Mechanics/Generalized Born Surface Area (MM-GBSA) approach was chosen for rescoring complexes since it is the quickest force field-based technique that estimates the free energy of binding, in comparison to the different computational free energy approaches, like free energy perturbation (FEP) or thermodynamic integration (TI) procedures. Based on comparative studies, the MM/GBSA methodology performed much better than the MM/PBSA (Molecular Mechanics/Poisson Boltzmann Surface Area) methodology [41].

ADMET & Lipinski's RULES
The discovery of a new molecule to be a potential drug requires an investigation of pharmacokinetic properties. Hence, the evaluation of ADMET properties is considered the best method to check and choose the target molecules via in silico methods using the online services available at Pkcsm.

Absorption
A review of the literature indicated that molecules or drugs with results (>30% abs) have more potency in crossing the intestinal barrier [33]. Consequently, the results show that all the compounds studied showed high values; the percentage of human oral absorption varied between 84.54 and 97.607%, indicating that these compounds were considered well absorbed.

Distribution
Distribution values predicted that the bioactive molecule of interest, Cynaropicrin has a low value of − 0.47 in Blood-Brain Barrier Permeability, which designates that don't enter the brain. Thus, undesirable effects will be reduced, proving their efficacy and drug-likeness behavior [34].

Metabolism
Computational metabolic behavior indicates the chemical biotransformation of a drug in the body; CYP plays a significant role in converting drug compounds. In Table 4, Cynaropicrin compound is revealed as substrates and inhibitors of the two main subtypes, 2D6 and 3A4 of   CYP, meaning they may have the ability to be metabolized in the liver [32][33][34].

Excretion and toxicity
The excretion result on Cynaropicrin shows an acceptable negative value; it is then evaluated for the Ames toxicity test, which is considered non-toxic; in parallel, in terms of Hepatotoxicity, Cynaropicrin has proven to be non-toxic, which justifies its effectiveness.

Lipinski's rule
The results of the toxicity tests of the ADMET method have proven that the bioactive molecule which interests us, Cynaropicrin, is the least toxic and the safest possible. Moreover, it can be observed from Table 5 that these three compounds comply with Lipinski's rule standards: they formed ≤5 H-bond donors and ≤10 H bond acceptors, MW < 500 Daltons, and octanol/water partition coefficient <5. The ADMET analysis and Lipinski's rules indicate that the designed Cynaropicrin was verified in silico as a safe pharmaceutical compound.
The compounds showed drug-like characteristics based on the rules of Lipinski and Veber (Table 5), which meant that they had a higher probability of good absorption after oral administration. According to Lipinski's rule, Veber and ghose these molecules exhibited drug-like properties.

Conclusion
Traditional therapeutic medicine is considered a real assistant for human health and well-being. In this research three bioactive molecules of Saussurea costus Costunolide, Cynaropicrin, and Dehydrocostus Lactone are studied and discovered as a potential medication for COVID-19. The reverse docking was conducted to determine the inhibitory action of these compounds' constituents with the different proteases of        Cynaropicrin and the examination of their in vitro activity against SARS-CoV-2 could be interesting and help to obtain new effective drugs.

Declaration of competing interest
1. We have no conflicts of interest to disclose.