In silico Screening of Potential SARS-CoV-2 Main Protease Inhibitors from Thymus schimperi

Abstract

Background

COVID-19 is still instigating significant social and economic chaos worldwide; however, there is no approved antiviral drug yet. Here, we used in silico analysis to screen potential SARS-CoV-2 main protease (M^pro) inhibitors extracted from the essential oil of Thymus schimperi which could contribute to the discovery of potent anti-SARS-CoV-2 phytochemicals.

Methods

The absorption, distribution, metabolism, excretion, and toxicity (ADMET) profiles of compounds were determined through SwissADME and ProToxII servers. AutoDock tools were used for molecular docking analysis studies, while Chimera, DS studio, and LigPlot were used for post-docking studies. Molecular dynamic simulations were performed for 200 ns under constant pressure.

Results

All compounds exhibited a bioavailability score of ≥0.55 entailing that at least 55% of the drugs can be absorbed unchanged. Only five (9%), nine (16%) and two (3.6%) of the compounds showed active hepatotoxicity, carcinogenicity, and immunotoxicity, respectively. Except for flourazophore P, which showed a little mutagenicity, all other compounds did not show mutagenic properties. On the other hand, only pinene beta was found to have a little cytotoxicity. Five compounds demonstrated effective binding to the catalytic dyad of the SARS-CoV-2 M^pro substrate binding pocket, while two of them (geranylisobutanoate and 3-octane) are found to be the best hits that formed hydrogen bonds with Glu¹⁶⁶ and Ser¹⁴⁴ of SARS-CoV-2 M^pro.

Conclusion

Based on our in silico analysis, top hits from Thymus schimperi may serve as potential anti-SARS-CoV-2 compounds. Further in vitro and in vivo studies are recommended to characterize these compounds for clinical applications.

Keywords:

• structural analysis
• SARS-CoV-2
• main protease
• inhibitors
• Thymus schimperi

Introduction

Coronavirus disease 19 (COVID-19) is caused by the novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Since its emergence in December 2019, it exerts a big economic and social impact globally. Data from worldometer (https://www.worldometers.info/coronavirus/) showed that COVID-19 cases surpass 655 million with a death toll of over 6.6 million as of December 14, 2022. It is unknown when the challenge posed by the pandemic will last. However, its impact could be longer due to various means of transmission, large asymptomatic carriers, absence of point-of-care cost-effective tests, absence of therapeutics and resistance of variants to vaccines.^1–6 One of the challenges in tackling the prevention and control activities is the presence of various modes of transmission.¹

Studies have been discovering promising inhibitors of SARS-CoV-2. Previously, the FDA approved remdesivir to use in case of emergency⁷; however, it has limited clinical outcomes in non-mechanically ventilated severely ill patients,⁸ indicating that further efforts are required to identify therapeutically effective antivirals.⁹ Owing to their key role in the viral cycle and relatively conserved nature, targeting viral enzymes like the main protease or 3C-like protease (M^pro or 3CL^pro), papain-like protease (PL^pro), non-structural protein 12 (nsp12), and RNA-dependent RNA polymerase (RdRP) could be promising.¹⁰

The main protease (M^pro) of SARS-CoV-2 plays key roles in viral replication and assembly^11,12 and human protease enzymes with similar specificity have not been reported so far,^13,14 making it the ultimate potential drug target. The M^pro is a cysteine protease where two protomers (each containing Domains I, II, and III) form a homodimer. Residues 8–101 and 102–184 form Domains I and II, respectively. A Cys-His catalytic dyad, which is proteolytically active, is located between Domains I and II.^15–19

Viral protease enzymes have a relatively conserved sequence, and they also have crucial roles in viral maturation and assembly which make them potential therapeutic targets. So far, the FDA approved several antiviral drugs targeting viral proteases including HIV-1 protease inhibitors²⁰ and hepatitis C virus (HCV) NS3/4A protease inhibitors.²¹ Since it is a key viral enzyme, the main protease of coronaviruses is a potential drug target that several studies are focusing on.¹¹ The M^pro digests polyproteins to release enzymes essential for replication and it also has NTPase and RNA helicase activity.^22,23 In this regard, studies looking for drugs against COVID-19 have been targeting SARS-CoV-2 M^pro. Therefore, it is essential to design antivirals targeting SARS-CoV-2 M^pro as it could have a potential clinical application.²⁴

Repurposed broad-spectrum drugs, new drugs, medicinal plants, and known antivirals effectively inhibited SARS-CoV-2 M^pro with potential antiviral effects.⁹ Among others, peptidomimetic alpha ketoamide inhibitors,^25,26 Michael acceptor compounds,²⁷ carmofur,^27,28 ebselen,^27,29 aldehyde-based compounds,³⁰ and 6e,³¹ clinically approved drugs (lopinavir/ritonavir),³² antiplatelet drug dipyridamole,^33,34 boceprevir (anti-Hepatitis C virus), GC-376, calpain inhibitors (II, XII) and GC-373^35–37 exhibited effective anti-SARS-CoV-2 M^pro activities.

Molecular docking and molecular dynamic simulations are becoming famous means of in silico drug discovery. These methods determine the effective binding of lead compounds in the active site of the M^pro of SARS-CoV-2. Accordingly, a study identified 12 best hits from Super Natural II and Traditional Chinese Medicine databases.³⁸ Besides, several potential compounds were identified from the traditional Chinese Medicine database interacting with active site residues of the M^pro (His⁴¹, Gly^143, and Cys¹⁴⁵).³⁹ Qumar et al⁴⁰ also recommended the top nine hits from several phytochemicals to consider as potential anti-SARS-CoV-2 M^pro lead molecules. Studies also showed that phytochemicals demonstrated strong pharmacokinetic and drug-likeness properties with acceptable toxicity profiles.⁴¹ Bioflavonoids⁴² also showed promising drug candidacy for SARS-CoV-2. In addition to targeting M^pro, phytochemicals have been demonstrating effective in silico drug potential for other core viral components including RNA polymerase.⁴³

The structure–function-based designing of drugs is an important phase in drug discovery. In addition to drug repurposing, studying the in silico structural basis of prospective SARS-CoV-2 M^pro inhibitors from medicinal plants is considered pivotal to further study the antiviral activity of these compounds for clinical application. In this regard, in silico mutational analysis of the M^pro helps for better drug repurposing.⁴⁴ Several studies have been investigating the in silico potential inhibitory activities of compounds against the M^pro of SARS-CoV-2. A study identified that three clinically approved drugs (glibenclamide, bedaquiline, and miconazole) effectively bind on the active site of SARS-CoV-2 M^pro with possible inhibitory activities.⁴⁵

In another study, Kanhed et al⁴⁶ identified ritonavir, nelfinavir, and saquinavir to be potent M^pro inhibitors. Ritonavir formed hydrogen bonds with Gly¹⁴³ and Cys¹⁴⁵ with its (thiazol-5-yl) methylcarbamate of oxygen, while the thiazolyl ring formed polar contacts with Thr²⁵, Thr²⁶ and Leu²⁷ of the S1’ subsite. Nelfinavir stabilized its binding with M^pro via hydrogen bonding with Glu¹⁶⁶, and with His⁴¹ and Tyr⁵⁴ in the S2 subsite. Depending on the structure of the pocket, compounds containing oxirane rings are suggested to be good M^pro inhibitors.⁴⁷ Arbutin, terbutaline, barnidipine, tipiracil, and aprepitant were identified as potential hits forming different hydrophilic, hydrophobic, and electrostatic interactions with M^pro.⁴⁸ Thioflavonol is a synthetic flavonoid analog that showed a strong binding with the conserved residues in the S1 subsite.⁴⁹ These and other studies illustrated the efficient binding of compounds on the substrate binding pocket (active) site of the M^pro indicating their potential anti-SARS-CoV-2 roles.

The genus Thymus contains about 350 species traditionally being used to treat different diseases worldwide.⁵⁰ Members of the genus, mostly T. schimperi, locally known as “Tosign”, in Ethiopia are used to treat diseases. Its leaves are used as spices in various Ethiopian foods and traditional medicines to treat different diseases including colds.⁵¹ No significant in vivo and silico rat embryonic toxicity was observed from compounds derived from the essential oil of this plant,⁵² indicating the further application of this plant extract for other infectious diseases including COVID-19. Here, we used in silico analysis to identify potential SARS-CoV-2 M^pro inhibitors from Thymus schimperi which could provide structural insights to discover potent anti-SARS-CoV-2 phytochemicals.

Materials and Methods

ADMET Analysis of Compounds

Extraction of 56 compounds from the essential oil of Thymus schimperi was done by using GC-MS.⁵² Then, the PubChem CID number of each compound was obtained from PubChem⁵³ and two-dimensional structures were built by Chemdraw (8.0).⁵⁴ SMILES were created through the Swiss ADME web tool.⁵⁵ Finally, the in silico absorption, distribution, metabolism, and excretion profiles were analyzed by using the online SwissADME server⁵⁵ and ProToxII predicted the in silico toxicity profile of the ligands.⁵⁶

Ligand and Protein Preparation

The ligands were prepared in PDB format first and then in PDBQT format suitable for molecular docking. Regarding protein preparation, the water molecules were removed, polar hydrogens were added, Kollam united atom force field was used to add charges and finally the structure was saved in a PDBQT format for further molecular docking analysis. To run the molecular docking, the crystal structure of SARS-CoV-2 M^pro (PDB ID: 6y2e)²⁵ was prepared and refined similar to previous studies by using AutoDock tools 1.5.7.^57–59 His⁴¹ was considered as the active site residue for molecular docking purposes. The protein we used (PDB ID: 6y2e) was solved by X-ray crystallography at a resolution of 1.75 Å.

Receptor Grid Preparation

The grid box for the M^pro was generated through AutoDock tools⁵⁸ around the active site. The Grid Box center was set as 75, 75, and 75 for the X, Y, and Z centers, respectively. The substrate-binding site is in the cleft between domains I and II and the protomers, which bind each other through N-terminus residues 1–7, are located between domains II and III with roles in the formation of the substrate-binding site.^{18,25,40,60,61} M^pro is a cysteine protease digesting viral polyproteins. First, a proton is transferred from Cys¹⁴⁵ to His⁴¹ with a simultaneous nucleophilic attack of the carbonyl carbon atom of the peptide bond by the sulfur atom of Cys¹⁴⁵ resulting in a thiohemiketal intermediate with subsequent protease activity.⁶²

Comparative Molecular Docking Analysis

Three docking tools (AutoDock Vina,⁵⁸ GOLD (Genetic Optimization for Ligand Docking) and MOE (Molecular Operating Environment)) were used for the molecular docking study of the 56 compounds (Supplementary Table 1). The selected compounds were sketched and minimized by ChemDraw Ultra 3D and saved into Mol2 and PDB format for molecular docking purposes.

The molecular docking of selected compounds with the M^pro was performed by AutoDock Vina to assess the inhibitory action of ligands against the M^pro. The ligand and protein were converted into pdbqt files using MGL tools. Grid sizes were adjusted to 75 × 75×75 Å in the X, Y and Z axes, respectively, with the grid spacing value of 0.650 Å to cover the target receptor. The fitness score was determined to select the best docking pose of inhibitors in the binding pocket. The best binding affinity was selected in all the docking processes.

After importing all the ligands into MOE, MDB file conversion, 3D protonation and energy minimization was conducted to 0.01 gradient. The protein structure concerns of the M^pro PDB structure were resolved in MOE. Before energy minimization, hydrogen atoms were introduced to structures with their typical geometry, and all solvent molecules were removed. Adjustments were made to the scoring methodologies’ default values. After the docking processes were completed, the acquired poses were analyzed, and the poses with the lowest allowable rmsd refine values and the same binding mode as the native ligand were chosen.

Alkyl, hydrogen, and hydrophobic interactions between the ligand and receptor were investigated through Chimera v1.8.1 and LigPlot⁶³ tools within the range of 5Å while BIOVIA Discovery Studio Visualizer v21.1 was used to identify the key interacting residues with maximum binding energies. Molecular dynamics (MD) simulations were conducted to validate docking results and assess the binding behavior and stability of prospective drugs by using the GROMACS 2016.4 Linux command line suite.²⁵

Validation of the Molecular Docking Method

The molecular docking protocol was validated by using an N3 peptide inhibitor as a reference ligand.²⁷ Then, a decoy set of ligands were used along with the active ligands. After the generation of the decoy dataset, all the ligands were prepared by MGL tool for docking.

Molecular Dynamics (MD) Simulation

The best-scored docking models of the most promising lead compound (3-octane) in complex with M^pro were selected as the starting coordinates for a 200-ns all-atom molecular dynamics simulation using the GROMACS-2016 program (GNU, General Public License; http://gromacs.org) and CHARMM36 force field. Solvation of the ligand–protein complex was done within a cubic box of the transferable intermolecular potential with a three-point (TIP3P) water model^93,93 with a minimum distance of 10 between the protein and each side of the 3D boxes. CHARMM General Force Field (CGenFF) tool was used to define the CHARMM force field parameters for the examined ligand⁶⁴ (ParamChem project; https://cgenff.umaryland.edu/). The protein residues were allocated to their standard ionization states under physiological conditions (pH 7.0), and the whole complexes were neutralized by adding enough K⁺ and Cl⁻ ions via the Monte Carlo ion-placing method. Heavy atom retaining, and maintaining original protein folding was also considered using a 1000 kJ/mol.N m² force constant in a three-step MD simulation. First, each system geometry was optimized with 5000 iterations (5 ps) of the steepest descent algorithm. In the succeeding step, the system was conditioned for 100,000 iterations (100 ps) at each stage of two-staged equilibration. The initial equilibration step was carried out using a constant number of particles, volume, and temperature (NVT) ensemble following the Berendsen temperature coupling method for controlling the temperature within the 3D box. Subsequently, the second equilibration stage was conducted in a constant number of particles, pressure, and temperature (NPT) ensemble at 1 atm and 303.15 K utilizing the Parrinello-Rahman barostat as a reference. Lastly, MD simulations were performed for 200 ns under constant pressure.

Results

ADME and Toxicity of the Compounds

Based on online server ADMET analysis, potential SARS-CoV-2 M^pro inhibitors from Thymus schimperi exhibited variable water solubility, GI absorption, blood–brain barrier (BBB) permeant, bioavailability score and synthetic accessibility. Among the screed compounds, 83.9% (47/56) of them are water-soluble, while the rest are either poorly soluble or moderately soluble. About 48% (27/56) of the compounds are predicted to be highly absorbable in the gastrointestinal (GI) system but 17.9% (10/56) of them are not able to penetrate the BBB. All compounds have a bioavailability score of ≥0.55 entailing that at least 55% of the drugs can be absorbed unchanged. Besides, almost all compounds have a synthetic accessibility (SA) score of below 5.00 while 35.7% (20/56) of them have a SA score of below 3.00 indicating their ease of synthetic accessibility (Supplementary Table 1).

The lethal dose (LD₅₀) of the compounds ranges from 113 to 5700 mg/kg. Two compounds, alpha-phellandrene and myrcenol, showed non-toxic (with a toxicity class value of 6 and LD₅₀ value >5000 mg/kg). Except for 6 (10.7%) compounds exhibiting a toxicity class of 3 indicating “toxic if swallowed”, other 48 (85.7%) compounds have a toxicity class of 4/5 indicating that these compounds are harmful/may be harmful if swallowed. Only five (9%), nine (16%) and two (3.6%) of the compounds exhibited active hepatotoxicity, carcinogenicity and immunotoxicity, respectively. Except for flourazophore P with a little mutagenicity, all other compounds did not show mutagenic properties. On the other hand, only pinene beta was found to have a little cytotoxicity (Supplementary Table 2).

Molecular Docking

In molecular docking, five compounds from Thymus schimperi (amphotericin-gamma, geranylisobutanoate, 3-octane, vetivenene_beta and germacrene-D) were the top hits based on lowest binding energy and highest gold score (Table 1). These compounds form hydrogen bonds, van der Waals interactions and alkyl interactions with the amino acids of M^pro (Supplementary Figure 1). Among these, amphotericin-gamma, geranylisobutanoate, 3-octane and germacrene-D showed strong binding on the His41-Cys-145 catalytic dyad of the SARS-CoV-2 M^pro substrate binding pocket (Supplementary Figure 2). Except for beta-vetivenene, other compounds exhibited strong interaction with the M^pro residues through all Chimera, DS studio and LigPlot analyses. Further, according to DS studio and LigPlot analysis, geranylisobutanoate, and 3-octane form hydrogen bonds Glu¹⁶⁶ and Ser¹⁴⁴ indicating that these compounds are the best among others (Table 1). Based on the highest binding affinity (−5.0 kJ/Mol) and gold score (44.88), 3-octane–M^pro protein–ligand complex was taken for MD simulation studies for investigating the potentiality of a small drug-like molecule to occupy the catalytic dyad of M^pro cavity in a way that would disrupt the function of M^pro (Figure 1).

Table 1 Summary of Top Five Hits Screened Against SARS-CoV-2 M^pro

Ligand Name	Binding Affinity (kcal/Mol)			Visualization
	AV	MOE	Gold Score	Chimera	DS Studio		LigPlot
				Interacting Residues	Interacting Residues	H Bond	Interacting Residues	H Bond
Geranylisobutanoate	−4.8	−5.17	44.26	His41, Met49, Phe140, Leu141, Asn142, Gly143, Ser144, Cys145, His163, His164, Met165, Glu166, His172, Gln189	His41, Met49, Phe140, Leu141, Asn142, Gly143, Ser144, Cys145, His163, Met165, Glu166, Gln189	1: Glu166	His41, Met49, Leu141, Asn142, Ser144, Cys145, His163, Met165, Glu166, Gln189	2: Glu166, Ser144
Amphotericin-gamma	−5.5	−4.8	37.98	Thr25, Thr26, Leu27, His41, Met49, Phe140, Leu141, Asn142, Gly143, Ser144, Cys145, His163, His164, Met165, Glu166, His172	Leu, 27, His41, Met49, Phe140, Leu141, Asn142, Gly143, Ser144, Cys145, His163, His164, Met165, Glu166		His41, Phe140, Leu141, Asn142, Gly143, Ser144, Cys145, His163, His164, Met165, Glu166
Germacrene-D	−5.0	−4.6	35.78	Thr25, His41, Ser46, Met49, Phe140, Leu141, Asn142, Gly143, Ser144, Cys145, His163, His164, Met165, Glu166, His172, Arg188, Gln189	His41, Met49, Asn142, Gly143, Cys145, His163, Met165, Glu166, Gln189		His41, Met49, Asn142, Cys145, His163, Met165, Glu166, Gln189
3-Octane	−5.0	−5.17	44.88	His41, Met49, Phe140, Leu141, Asn142, Gly143, Ser144, Cys145, His163, His164, Met165, Glu166, His172, Gln189	His41, Met49, Phe140, Leu141, Asn142, Gly143, Ser144, Cys145, His163, His164, Met165, Glu166, His172, Gln189	1: Glu166	His41, Met49, Leu141, Asn142, Ser144, Cys145, His163, Met165, Glu166, Gln189	2: Glu166, Ser144
Vetivenene-beta	−5.6	−4.9	36.88	Glu55, Leu58, Ile59, Lys61, His80, Ser81, Met82	Glu55, Leu58, Ile59, His80, Ser81, Met82

Figure 1 3D and 2D representation of the best docked pose of the compound 3-octane highlighting the critical binding site residues of M^pro. (A) Ribbon representation of the ligand 3-octane bound to the M^pro. The blue color represents the ribbon representation of M^pro while the ligand is represented in green color. The protein residues interacting with the ligand are represented in orange color. (B) The surface representations of the M^pro depicting the ligand are tightly bound within the binding groove of the protein. (C) The 2D representation of the docked pose protein–ligand complex. (i) 2D interaction diagram of 3-octane–M^pro complex generated by MOE. (ii) Ligplot representation of the protein–ligand complex drawn by LigPlot+. (iii) 2D depiction of the docked complex highlighting hydrogen bonds and cation-p interaction through DS Visualizer.

Molecular Dynamics Simulation Studies

The stability of the projected docked 3-octane–M^pro complex was evaluated by an all-atom MD simulation study. Implementing this type of study would also give useful insights concerning the dynamic nature of both the ligand and M^pro, as well as an evaluation of the ligand’s most significant binding interactions with essential catalytic site residues.⁶⁵ Enrolling of the predicted ligand–protein complex for 3-octane and M^pro protein was done within a 200 ns MD simulation run.

Trajectory Analysis of 3-Octane–M^pro Complex

MD experiments conducted in the complex elucidated the plausible mechanism of inhibition. The shape of the protein influences the conformational dynamics, thus understanding the functional flexibility of a biological macromolecule is of utmost important.⁶⁵ The average RMSD for the 3-octane–M^pro complex was found to be 0.25 Å when computed for 200 ns which indicated that the system was highly stable during the MD simulations (Figure 2A). Similar to previous RMSD analysis, the obtained Rg plot confirmed better protein–ligand complex stability with an average value of 0.25 Å (Figure 2B).

Figure 2 Global stability analysis of ligand-hACE2 protein complexes throughout 200 ns all-atom MD simulation. (A) Analysis of RMSD trajectories for the 3-octane–M^pro protein complex throughout 200 ns MD simulation deciphering the primary conformational switches. (B) The radius of gyration for the protein–ligand complex reflects the complex structure’s global stability. (C) RMSF of protein–ligand complex depicting fluctuations across protein residues during 200 ns of MD simulations. (D) Extent of M^pro binding site coverage via SASA analysis along with the time evolution 200 ns all-atom MD simulation.

As shown in Figure 3C, RMSF values were 0.2 Å indicating the flexibility of residues forming the catalytic dyad. The 200 ns of MD run showed that the flexibility of the residues reached a peak upon binding of the ligand with the receptor as depicted by the RMSF. Significant fluctuations were observed at SER-46, LEU-50, THR-190, THR-224, and ASN-277 and PHE-305 with an average RMSF value of 0.2 Å (Figure 2C). SASA corresponds with the molecular surface area that can be assessed by solvent molecules, hence providing a quantitative assessment of the level of protein/solvent interaction. The study was performed on the atoms of residues bordering the M^pro binding site to estimate the solvent-exposed region (Figure 2D).

Figure 3 The structure superimposition of MD trajectories (0 ns, 50 ns, 100 ns, 150 ns, and 200 ns) over the pre-simulated 3-octane–M^pro complex. (A) The secondary structure representation of superimposed protein–ligand complexes. The ligand is represented in green color. (B) The surface representation of M^pro representing that ligand is bound within the same active site groove throughout the 200 ns of MD simulation run considering it stable active site conformation of the ligand.

Conformational Analysis Across the Selected Trajectories

To detect the structural variability of the simulated system, the MD trajectories (0 ns, 50 ns, 100 ns, 150 ns and 200 ns) were superimposed on the pre-simulated 3-octane–M^pro docked complex (Figure 3) and produced the superimposed structure at 0.8 Å. The 3D and 2D interaction diagrams (Figure 4) of 0 ns, 50 ns, 100 ns, 150 ns and 200 ns trajectories presented that the interacting residues of M^pro were consistent with the course of the 200 ns MD simulation run. Interestingly, there is no significant orientation change for the ligand within the M^pro binding site between the time frames 0 and 200 ns. A slight shift caused a loss of the initial hydrogen bond with GLU-166. Stabilization of 3-octane within its new conformation/orientation was further mediated by several hydrophobic residues including His⁴¹, Leu¹⁴¹, Ser¹⁴⁴, and Cys¹⁴⁵.

Figure 4 Conformations of the 3-octane complex at M^pro binding site through selected trajectories.

Discussion

Computational drug design is quite an important tool to speed up drug discovery, especially in case of emergency.⁶⁶ So far studies focused on both repurposing previously approved drugs and/or identifying new potent phytochemicals. Proposed plant metabolites are assumed to serve as potential anti-SARS-CoV-2 lead compounds to combat COVID-19.⁶⁷ Medicinal plants especially those used in Chinese traditional medicine have been investigated for applications in COVID-19 treatment.⁶⁸ Studies identified potential SARS-CoV-2 M^pro inhibitors from different plants and Traditional Chinese Medicine databases that effectively interact with the active site residues in the catalytic dyad of the Mpro^{38–40,69–73} with promising outcomes.

Here we screened 56 compounds from Thymus schimperi for their ADME, toxicity and binding profiles. Most of the screened compounds showed good water solubility, gastrointestinal absorption, and permeant to the BBB. According to bioavailability score definitions,⁷⁴ our compounds could be well absorbed in humans since they have a good predicted bioavailability score of ≥0.55 which is similar to previous studies.⁷⁵ Synthesizing natural compounds could be better to minimize the resource and time required to extract from their natural sources. In this regard, our molecules have an SA score below 5.00 indicating that these compounds are not difficult to synthesize,⁷⁶ indicating their potential consideration for further optimization and drug development studies.

Five compounds from Thymus schimperi (amphotericin-gamma, geranylisobutanoate, 3-octane, vetivenene-beta and germacrene-D) were the top hits that formed hydrogen bonds, van der Waals interactions and alkyl interactions with the amino acids of M^pro. Specifically, geranylisobutanoate and 3-octane showed strong SARS-CoV-2 M^pro binding affinity with acceptable toxicity and ADME profiles except these compounds showed hepatotoxicity in in silico analysis. Interestingly, these compounds are under toxicity class 5 which is acceptable in medicine. The effective binding of these compounds to the M^pro and their good ADMET profiles is in line with previous studies.^77–82 Thus, considering these compounds for further applications is of significance; however, in vitro, and in vivo chronic toxicity studies are pivotal.⁵²

Previous studies also reported similar findings to our study. Numerous plant-derived phytochemicals such as curcumin, gartanin, robinetin,⁸³ amentoflavone, gallocatechin gallate,⁸⁴ chelidimerine, rutin, fumariline, catechin gallate, adlumidine, astragalin, somniferine,⁸⁵ kaempferol, herbacetin, eugenol, 6-shogaol,⁸⁶ triacontane, hexacosane, methyl linoleate, and methyl palmitoleate,⁸⁷ cosmosiine, pelargoniding-3-O-glucose, and cleomiscosin⁸⁸ formed similar binding patterns while docking with the main protease. Another study also reported that plant-derived phytochemicals including flavan 3-ols (catechins/procyanidins), complex oligomeric procyanidins (procyanidin A3, procyanidin A4, procyanidin A1, and procyanidin B3) exhibited very good binding characteristics like our compounds.⁸⁹ Interestingly, the main protease of coronaviruses is relatively conserved.^90–92 Thus, designing drugs against the main protease could broadly attack the virus variants. Additionally, the absence of human host-cell proteases with similar activity with viral main proteases make the prospective drugs to have less off-target actions.^13,14

Bioinformatics analyses have demonstrated their value in the invention of innovative computer-assisted compounds against a variety of diseases, including neurological disorders, cancer, and pathogenic infections.⁹³ To better comprehend the structure and function of M^pro, a crucial element in the design of medications, computational studies are necessary. Computational drug design methods also play a crucial part in determining which medicine is the best among others.²⁴ MD simulation studies connect the theoretical and experimental studies in drug discovery by analyzing the protein–ligand complex. Further, these studies are crucial to determine the stability of the ligand-receptor complex.⁹⁴ The MD simulation results in this study confirm the firm binding of 3-octane at 200 ns. These results present pieces of evidence that 3-octane could be a promising SARS-CoV-2 M^pro inhibitor, while in vitro and in vivo studies are yet required for further clinical application. The problem behind drug discovery is less probability of the anticipated compounds passing clinical trials. Therefore, continuous efforts and further investigation of compounds through in silico, in vitro and in vivo studies are required.

Compounds from food spices are reported to have potential anti-SARS-CoV-2 M^pro activities.⁹⁵ Thymus schimperi is used as a food spice in Ethiopia. Studies have shown that essential oils from Thymus schimperi exhibited promising antimicrobial activities.^96,97 In silico screening of lead compounds from food promote the therapeutic role of food in combating the COVID-19 pandemic.⁹⁸ As the top hits in this study exhibited excellent binding and acceptable ADMET properties, further studies are recommended to optimize these compounds in the plant to design drugs tackling respiratory diseases. More importantly, considering geranylisobutanoate and 3-octane for future clinical applications could be of significance.

Many lessons have been learned from previous in silico and structure-based drug designing studies that would help prospective studies succeed fast in discovering effective antivirals to tackle COVID-19 and other related diseases.²⁴ Future studies need to focus on the inhibitor enzyme complex which included atomic-level mechanisms of peptide cleavage, pharmacophore requirements of the M^pro, stability of the inhibitor-enzyme complex, and plasticity of the active site of M^pro. Besides, the occurrence of mutations at the domains and/or the active site affecting the pocket, the size and accommodation capacity of the subsites of the M^pro should be considered in designing new drugs or modifying previously known broad-spectrum drugs. Most studies solely report the binding affinity and energy of compounds towards the substrate-binding cleft of the M^pro, improvements considering the abovementioned points should be taken in the future.

Generally, structure–function analysis is a key factor in drug design and discovery where in silico studies are crucial to predict the drug ability of compounds in a short period in addition to its indispensable role in predicting the best drug candidates among others. The big challenge in computational drug designing is that the clinical use of these desired drugs is questionable corresponding to the possible limitations of passing clinical trials. Although our study characterizes a limited number of compounds, it provides a glimpse into the identification of lead compounds from food spices that could serve as therapeutic drugs in clinical medicine.

Conclusion

We identified potential compounds from Thymus schimperi that could serve as drugs to inhibit SARS-CoV-2 M^pro. Among the top hits, geranylisobutanoate and 3-octane exhibited strong binding to the catalytic dyad of SARS-CoV-2 M^pro with acceptable ADMET profiles. Interestingly, our results entail that 3-octane could be the most promising anti-SARS-CoV-2 M^pro inhibitor. Further studies are required to determine the in vitro and in vivo characteristics of the identified potential compounds to employ them in clinical applications.

Ethical Approval

Our study is a pure computational study and does not include human data and thus no need of ethical approval.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising, or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Disclosure

The authors declare that there are no conflicts of interest.