Phytochemicals

A total of 129 phytocompounds from 34 Algerian aromatic and medicinal plants were selected. The plants were chosen on the basis of their traditional medicinal use and the amount and type of active principles present in them. The first plants selected were those used to treat respiratory infections and illnesses.

Software programmes

In this study, ChemSpider and Drugbank were used as the databases to find all molecular structures ^[28-29]. The RCSB Protein Data Bank's crystal structure for the receptor was used (PDB code 6LU7) ^[30]. All the preparation was done using Autodock tools 1.5.6 ^[31]. Autodock vina was used to carry out the docking study ^[32,33]. The pkCSM ^[34] and Swissadmet ^[35] online servers were carried out to calculate the Lipinski's rule of five and ADMET properties ^[36]. Simulation of molecular dynamics MD was carried out using Gromacs 2019 package ^[37] and CHARMM 2017 force field. The results were analysed applying the Discovery studio 2016 package ^[38].

Molecular docking

Molecular docking, which is performed prior to conducting any experiments, can show that any biochemical process is feasible ^[39]. Generally, there is four-step procedure to perform the molecular docking, namely, target selection, Ligand selection, preparations, and docking. In this study, molecular docking was performed to predict poses in which the different phytocompounds may be in the Mpro active site. This allows analysing possible intermolecular interactions between the target and the ligand, where various orientations and conformations of the protein interacting with its ligand are generated ^[20].

The extended PDB format (PDBQT) was employed for ligand and enzyme preparations. Further receptor preparation includes the elimination of cofactors, water molecules, and other ligands; and the addition of polar hydrogens. To obtain the binding site, the crystallographic structure of SARS-CoV-2 Mpro (PDB code 6LU7) was imported into the programme's "work space" in Discovery Studio 2016 ^[40]. Based on the co-crystallised inhibitor, the active site's centre was then identified and its coordinated extracted (x= -10.782, y= 15.787, and z= 71.277) ^[22]. For the x, y, and z points, the grid size was set to 20x20x20, and a grid spacing of 1 Á was established to cover the full binding pocket.

Lipinski's rule of five and ADMET prediction

The rule of five (ROF) is a general guideline for assessing drug similarity or determining whether a chemical compound with a particular pharmacological or biological activity possesses characteristics that would make it likely to be an orally active drug in people ^[41]. Molecular weight (≤ 500), logP (< 5), the number of rotatable bonds (< 10), the number of hydrogen bonds (HB) acceptor (≤ 10), and the number of HB donors (< 5) are among the parameters computed (values of a given property which will not be exceeded) for the Lipinski's rule of five; so, for more than one violation of the ROF, molecules will be viewed as being unimportant for future development, according the authors themselves.

In the study, ADMET was used to predict the absorption, distribution, metabolism, excretion, and toxicity parameters of the examined compounds.

Molecular dynamic simulations

Macromolecular structure-to-function relationships can now be effectively understood using a mature method called molecular dynamics simulations. Current simulation times are close to biological times and the wealth of data collected on the dynamic properties of macromolecules allows structural bioinformatics to change its standard paradigm from studying single structures to analysing conformational ensembles ^[42].

For this study, simulations of molecular dynamics (MD) were carried out only on the two best molecules derived from the docking studies. The main objective of MD simulations is to evaluate the stability of the ligand-enzyme complex found in the docking study. The results obtained from the docking calculation were used as input for the MD simulation. The ligand-receptor complex was prepared by adding water molecules to solvate the system and sodium and chlorine atoms to neutralise it. Prior to the simulation, the system was relaxed for 100 ps at 300 K and 1 bar. Then, simulation was run for 20 ns.

Molecular docking

Among the five therapeutic protein targets of SARS-CoV-2: 3CLpro (main protease), PLpro (papain-like protease), SGp-RBD (spike glycoprotein-receptor binding domain), RdRp (RNA dependent RNA polymerase), and ACE2 (angiotensin-converting enzyme 2) the poses and potential types of interactions between the 129 compounds and SARS-CoV-2 Mpro were carried out using molecular docking. The results are reported on Table 1. Among the 129 molecules docked, warfarin (compound No. 107) was previously studied by us ^[43].

Table 2 Structures that bind to SARS-CoV-2 Mpro with the greatest affinity. 

Warfarin has the highest docking score out of all the compounds (-7.1 kcal/mol). Here it is interesting to mention that warfarin is an anticoagulant and is used to prevent thrombosis and embolism. Its anticoagulant effects are based on blocking the vitamin K-dependent (VKD) stage in the full synthesis of several blood clotting components necessary for healthy blood coagulation ^[53]. In consequence, interactions between warfarin and Chinese herbal medicine were reviewed and described ^[54]. Additionally, some features of warfarin's anti-tumor effect are thought to be a result of its antithrombotic activity. It's also crucial to note that warfarin's inflammatory and immunomodulatory actions point to a far wider range of possible consequences ^[55]. Additionally, warfarin toxicity and metabolism are highlighted ^[56].

Other than warfarin, six additional molecules were found to have interesting docking scores which include catalponol, benzyl salicylate, Prasterone, davana ether, chamazulene, and rosifoliol (Table 2). Therefore, these seven molecules underwent further study for being the best possible candidates to inhibit SARS-CoV-2 Mpro.

Some of these compounds were found to possess biological activities, such as chamazulene, which is considered a natural anti-inflammatory ^[57] associated to its antioxidant properties, and contributes to the in vivo antiphlogistic activity by inhibiting the leukotriene synthesis, and it has antitumor activity as well ^[58].

According to Lobo ^[59], Prasterone has been approved by the US Food and Drug Administration (FDA) for symptoms of dyspareunia in postmenopausal women; and Heo ^[60] proposes intravaginal Prasterone as a treatment for vulvovaginal atrophy in postmenopausal women. Huang et al. ^[61] confirmed the pharmacological functions of catalponol which improved the biosynthesis of dopamine through the activation of tyrosine hydroxylase and the protection of PC12 cells cytotoxicity induced by L-DOPA.

To verify the reliability of the procedure and the software used in this work to produce reliable results, the co-crystallised ligand (N3) was re-docked into the active site and a docking score of -6.9 kcal/mol was found. Successively, the N3 structure was compared to the docked one (Figure 1a) and 16 interactions were found (Figure 1b ^[62]). The interactions between warfarin and the Mpro were established previously ^[43]. For the other interactions between compounds and SARS-CoV-2 main protease, they are shown in Figures 2 (a-f).

Figure 1 (a) N3 co-crystallised ligand (red) and N3 docked ligand (cyan); (b) Interactions between inhibitor N3 and SARS-CoV-2 Mpro^[53]. 

Figure 2 Interactions between Catalponol (a), Benzyl salicylate (b), Prasterone (c), Davana (d) Chamazulene (e), Rosifoliol (f) and SARS-CoV-2 main protease 

The results show that both structures overlap in most of the molecule, demonstrating that the procedure employed is able to reproduce the experimental results. The difference between the two structures is noted at one end of the molecule, where the docking result shows that the chain is oriented upward due to the flexibility of the peptide losing the interaction with the receptor, while the co-crystallised structure is interacting with the whole pocket. Analysing the N3-Mpro complex, 16 interactions were found, which include a Van der Waals interaction with Asn142, HB interactions with Glu 166, His 164, Gly 143, Thr 190, Gln 189, His 163, and Phe 140 residues, amide-n Stacked interactions with Leu 141 residue, Carbon HB with Met 165 and His 172 residues, and Alkyl and/or n-Alkyl interactions with Leu 167, Ala 191, Met 167, Pro 168, Met 49 and His 41 residues ^[62].

This result is relevant to analyse whether the studied molecules may interact with the active site and whether those interactions may be strong enough to compete with the natural substrate inhibiting the enzyme. Based on the docking score, only warfarin shows a higher value compared to the N3 inhibitor. The docking results of the best seven molecules based on the docking score show that Warfarin interacts with five residues ^[43], catalponol with five residues, benzyl salicylate with six, prasterone with two, davana ether with three, chamazulene with three, and rosifoliol with two residues. It is important to notice that the N3 inhibitor forms more than twice the interaction of any of these drug candidates which may infer a weak interaction in case any interaction at all is produced. Delving deeper into the analysis of the interactions between warfarin and the Mpro, the interactions include HB with Glu 166, Ser 144 and Cys 145 residues, n-n T-shaped with His 41 residue and n-alkyl with Met 164 residue ^[43]. Also, the presence of HB and n-n T-shaped interactions in the complex formed by warfarin and Mpro increase the affinity of the complex and may be related with its pharmacological activity. This is in accord with the reports in the literature where the HB role is identified with the pharmacological effect of a ligand ^[21,62].

The values obtained for Mpro docking with warfarin, catalponol, benzyl salicylate, prasterone, davana ether, chamazulene and rosifoliol may be related to the number and type of bonds observed in the ligand-target formed, in fact the presence of the HB suggests an important potential pharmacological effect for the Mpro inhibition of SARS-CoV-2. Thus, the inhibition of this protein will stop viral replication and may become useful in the clinical management of COVID-19. These results are supported by many other studies where Glycyrrhiza glabra is used in traditional medicine to treat cough and as a laxative as mentioned by Chouitah ^[63]. Warfarin and catalponol, found in this study to be two of the most interesting molecules, are part of the active ingredients present in Glycyrrhiza glabra. Therefore, the extract of this plant could have very interesting effects on the inhibition of SARS-CoV-2 Mpro. According to Soleiman-Meigooni, liquorice (Glycyrrhiza glabra L.) syrup can be used as an adjuvant treatment for mild to moderate Covid-19 in addition to the recommended course of care ^[64]. It was found that, glyasperin A, glycyrrhizic acid, and compounds having oxane ring and chromenone ring substituted with hydroxyl 3-methylbut-2-enyl group could be considered as the best molecules and leads for the development of a new drug based on liquorice against COVID-19 ^[65,66].

Lipinski's Rule of five and ADMET Prediction

Using the web servers pkCSM and Swissadmet, the Lipinski's rule of five, and the ADMET parameters of the chemicals examined in this paper were computed. According to "the rule of 5", an orally active drug may have no more than one violation of the following criteria (Table 3): molecular weight less than 500; not more than 5 H-bond donors, 10 H-bond acceptors and 10 rotatable bonds; an octanol-water partition coefficient log P not greater than 5.

Table 3 Lipinski's rule of potential inhibitors. 

The results obtained for all the studied compounds in this work respect Lipinski's rule of five, including warfarin ^[43]. The N3 inhibitor has three Lipinski violations ^[30], which is to be expected, since it is a peptide-type inhibitor, characterised by a high MW, and many more than 10 rotational bonds. The parameters predicted (Table 4) were water solubility (WS) and human intestinal absorption (IA), volume of distribution (Vd), blood-brain barrier permeability (Bbp), cytochrome P450 (CYP450) metabolism, total clearance (TCl), hERG I and II inhibition (hERGi) and syntetic accessibility (SA).

Table 4 Computer simulation of ADMET forecast and accessibility values. 

All the molecules exhibit a low value for the blood-brain partition coefficient. This may indicate that, there will not be any risk of central nervous system (CNS)-related toxicity because these molecules have poor potential to pass the blood-brain barrier. The WS, IA, and metabolism pharmacokinetic parameters present drug-like values with more than 92% of intestinal absorption, WS between 6.3x10^-4 and 4.0x10^-6 mol/l, and a hepatic metabolism for catalponol, benzyl salicylate, and Prasterone. Also, the toxicity predicted for hERG I and hERG II inhibitors indicates that any compound exhibits hERG I, while only prasterone has hERG II inhibition. SA values show the feasibility of experimentally synthesising all the studied compounds.

Molecular dynamics simulation

Although the seven studied compounds were found to have good docking scores against Mpro, following Lipinski's rule of five and exhibiting appropriate ADMET properties; only the two best compounds (warfarin and catalponol) were taken to further MD simulations to determine the possibility of keeping the ligand-receptor complex in a solvated system for a period of time. Results will increase the feasibility of these compounds for being possible COVID-19 treatments through Mpro inhibition. In Figure 3, the results for the simulation of Warfarin against Mpro is shown.

Figure 3 (a) Potential energy; (b) RMSD of Warfarin (red) and the Warfarin-Mpro complex (gray); (c) number of HBs; (d) number of pairs within 0.35 nm; (e) Coulomb (blue) and Lennard-Jones (orange) interaction energy during the simulation. 

For this complex, the potential energy (Figure 3a) reaches a plateau before 500 ps showing that after this time, the system is at its lowest potential energy. The RMSD plot (Figure 3b) shows the ligand-enzyme complex maintains values lower than 0.5 nm while the ligand presents values lower than 1 nm during most part of the simulation (first 15 ns). During the last 5 ns the RMSD values in the ligands increase to around 1.1 nm. Looking at the HB formation and the number of pairs within 0.35 nm (Figures 3c and 3d) it may be observed that the maximum number of HBs and pairs is four and the minimum number is zero. In most parts of the simulation, at least one HB and one pair are formed. Finally, the Coulomb and Lennard-Jones interaction energies were evaluated (Figure 3e). While the Coulomb energies stayed constant throughout the simulation with values around -5 kcal/mol, the Lennard-Jones interaction energies increased as the simulation proceeded starting at -25 kcal/mol and ending at around -17 kcal/mol. Comparing the structure at the beginning of the simulation and at the end, it was observed that warfarin moves out of the active site cavity and stayed in the most superficial part of the pocket (Figure 4). This movement will make warfarin lose most of its potency although there is still some chance based on the interaction formed to be strong enough to compete with natural substrates.

Figure 4 Docking (Purple) and MD (Green) results of the interaction between warfarin and Mpro. N3 inhibitor is shown in pink. 

Analysing the simulation between catalponol and Mpro, the potential energy plateau is reached at 1181 ps. The three-dimensional structure of the complex appeared after 20 ns, the ligand could not manage to stay in the active site and leaves the system. This can be reflected in the plots found in Figure 5. The RMSD plot shows catalponol RMSD values are less than 0.5 nm in the first 5 ns of simulation. In the next 5 ns, values rise up to 2 nm. After that, in the final 10 ns, values experience a great increase reaching values of 8 nm. Looking at the HB and the number of pairs within 0.35 nm, catalponol makes interactions only during the first 10 ns of simulation. After that, all interactions are lost. In the first 10 ns, catalponol makes between 1 and 4 HBs and between 2 and 7 pairs within 0.35 nm. The same can be observed in the Coulomb and Lennard-Jones interaction energy where negative values are found up to 12 ns. After that all interaction is lost and both energies are zero.

Figure 5 (a) Potential energy; (b) RMSD of catalponol (red) and the catalponol-Mpro complex (gray); (c) number of HBs; (d) number of pairs within 0.35 nm; (e)Coulomb (blue) and Lennard-Jones (orange) interaction energy during the simulation.