Phytochemicals from Leucas zeylanica Targeting Main Protease of SARS-CoV-2: Chemical Profiles, Molecular Docking, and Molecular Dynamics Simulations

Simple Summary Molecular docking in conjunction with molecular dynamics simulation was accomplished as they extend an ample opportunity to screen plausible inhibitors of the main protease from Leucas zeylanica. The preferential phytochemicals were identified from L. zeylanica through gas chromatography–mass spectrometry (GC-MS). The pre-eminent three identified phytochemicals exhibited toxicity by no means during the scrutinization of ADME/T prominences. Moreover, pharmacologically distinguishing characteristics and the biological activity of the lead phytochemicals were satisfying as an antiviral drug contender. Additionally, the molecular dynamics simulation exhibited thermal stability and a stable binding affinity of the protein–compound complex that referred to the appreciable efficacy of lead optimization. Therefore, the preferable phytochemicals are worth substantial evaluation in the biological laboratory to recommend plausible antiviral drug candidates. Abstract Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a contemporary coronavirus, has impacted global economic activity and has a high transmission rate. As a result of the virus’s severe medical effects, developing effective vaccinations is vital. Plant-derived metabolites have been discovered as potential SARS-CoV-2 inhibitors. The SARS-CoV-2 main protease (Mpro) is a target for therapeutic research because of its highly conserved protein sequence. Gas chromatography–mass spectrometry (GC-MS) and molecular docking were used to screen 34 compounds identified from Leucas zeylanica for potential inhibitory activity against the SARS-CoV-2 Mpro. In addition, prime molecular mechanics–generalized Born surface area (MM-GBSA) was used to screen the compound dataset using a molecular dynamics simulation. From molecular docking analysis, 26 compounds were capable of interaction with the SARS-CoV-2 Mpro, while three compounds, namely 11-oxa-dispiro[4.0.4.1]undecan-1-ol (−5.755 kcal/mol), azetidin-2-one 3,3-dimethyl-4-(1-aminoethyl) (−5.39 kcal/mol), and lorazepam, 2TMS derivative (−5.246 kcal/mol), exhibited the highest docking scores. These three ligands were assessed by MM-GBSA, which revealed that they bind with the necessary Mpro amino acids in the catalytic groove to cause protein inhibition, including Ser144, Cys145, and His41. The molecular dynamics simulation confirmed the complex rigidity and stability of the docked ligand–Mpro complexes based on the analysis of mean radical variations, root-mean-square fluctuations, solvent-accessible surface area, radius of gyration, and hydrogen bond formation. The study of the postmolecular dynamics confirmation also confirmed that lorazepam, 11-oxa-dispiro[4.0.4.1]undecan-1-ol, and azetidin-2-one-3, 3-dimethyl-4-(1-aminoethyl) interact with similar Mpro binding pockets. The results of our computerized drug design approach may assist in the fight against SARS-CoV-2.


Introduction
In recent decades, the world has witnessed an unprecedented number of life-threatening human disease outbreaks caused by an array of pathogenic organisms, including several notable viral diseases, such as influenza, chikungunya, Nipah, Zika, and Ebola [1,2]. However, the ongoing spread of coronavirus disease 2019 (COVID-19) has been exponential. It has already surpassed most previous viral infections in terms of infectivity and has become the center of global attention. Wuhan, a populous Chinese city located in the Hubei province was the first location where this acute respiratory infection was identified in late December 2019 [3]. COVID-19 has taken a significant toll on people worldwide and on 11 March 2020 was declared by the World Health Organization (WHO) a pandemic. This highly contagious infection has had a detrimental impact on the global healthcare management system and, as of 1 February 2021, >100 million confirmed cases have been reported, including more than 2 million estimated deaths worldwide [4].
SARS-CoV-2 is a pleomorphic, enveloped, nonsegmented, single-stranded RNA betacoronavirus belonging to the Coronaviridae family and features a large genome (27-32 kb) that encodes both structural and nonstructural proteins. SARS-CoV-2 is associated with a higher transmission rate than other well-known human beta-coronaviruses, such as SARS-CoV and Middle Eastern respiratory syndrome coronavirus (MERS-CoV) [5]. The coronavirus main protease (M pro ) is a nonstructural protein that plays a crucial role in protein translation, viral replication, and maturation [6,7]. In a recent study, Liu et al. confirmed the existence of the M pro (also known as 3CL Pro or chymotrypsin-like protease) enzyme in SARS-CoV-2 [8]. The genome of SARS-CoV-2 encodes pp1a and pp1ab, two large polyproteins, similar to other Coronaviridae genomes [9]. The resulting polyproteins, pp1a and pp1ab, must be cleaved to generate mature nonstructural proteins (nsps) [10]. The large pp1a (replicase 1ab) is generated inside the cell via genomic RNA transcription. Therefore, the inhibition of M pro activity is anticipated to result in the prevention of viral replication, as similar cleavage specificity has not been identified in any human proteases [11], indicating that the polypeptide cannot be properly cleaved in the absence of M pro . Additionally, M pro has very low cytotoxicity and low similarity with human proteases [12]. Proteins produced by this pathogenic organism have been demonstrated to intervene the host immune response, and M pro enzyme-specific T cells have been encountered in SARS-CoV-2 patients [13,14]. Therefore, M pro is considered to represent a promising drug target for antiviral drug development.
In general, cytokine production, cell death, inflammation, and other pathophysiological processes are commonly associated with disruptions in redox balance, resulting in oxidative stress during viral infections, which can negatively affect the respiratory tract. Previously, viral replication was strongly correlated with the excessive production of reactive oxygen species (ROS) and a reduction in the components of antioxidant mecha-nisms [15][16][17][18]. Inflammatory reactions are triggered by the COVID-19 infection, resulting in the subsequent release of proinflammatory cytokines, which can cause acute lung damage [19]. Oxidative stress also plays a critical role in the perpetuation of the cytokine storm cycle and is important for blood clotting mechanisms [20]. The observed increase in COVID-19 infection severity in patients diagnosed with chronic diseases has been linked with the poor performance of the antioxidant system, suggesting that antioxidants may represent a prospective therapeutic option for COVID-19 infection [21].
A close connection exists between the innate immune response and the thrombotic response, and recent COVID-19 clinical data have revealed a correlation between this infection and thrombotic complications, which might result in increased incidence of microvascular thrombosis, venous thromboembolic illness, and stroke. Markers of COVID-19 include thrombotic complications, which are often associated with multiorgan failure and increased fatality [22].
Recently, phytochemicals have been investigated against different target proteins of SARS-CoV-2 to find appropriate lead compounds for COVID-19 infection. Baicalin, baicalein, 25-hydroxycholesterol, chrysosplenetin, shikonin, panduratin A, and quercetin are some of the examples of plant compounds that exhibited a potential effect against SARS-CoV-2 during in vitro studies. In silico studies were conducted on a broad spectrum for a plethora of medicinal plants and a huge number of compounds were screened. This not only helped indicate the potentiality of natural plant compounds but also reduced the number of tedious and costly wet-laboratory experiments [23]. Therefore, in our current study, we endeavored to assess the roles of phytoconstituents identified from Leucas zeylanica, a medicinal plant belonging to the Lamiaceae family, in the management of COVID-19 infection by employing computational biology approaches. Phytocompounds possess a wide range of pharmacological activities, and traditional healers have employed plants belonging to the Leucas genus to treat various disease states, indicating an immense potential for the discovery of new lead compounds [24]. L. zeylanica is a weed commonly referred to as "Ceylon slitwort" and locally known as "Kusha" [25]. The plant is widely distributed throughout China, India, Bangladesh, Sri Lanka, Thailand, Indonesia, Philippines, Vietnam, Cambodia, Nepal, Myanmar, Malaysia, and New Guinea [26]. A phytochemical screening of an L. zeylanica methanol extract confirmed the presence of alkaloids, flavonoids, tannins, steroids, and glycosides, which contribute to the traditional medicinal properties of the plant [27]. This plant is traditionally used as a vermifuge ingredient in addition to the treatment of burning sensations during urination, scabies, convulsion, fever, jaundice, scorpion and snake bites, colds, rheumatism, roundworm, psoriasis, anorexia, flatulence, colic, and malaria [24,25,28,29]. The antimalarial drugs chloroquine and hydroxychloroquine have been suggested as potential anti-COVID-19 therapies in recent studies; therefore, the traditional use of L. zeylanica against malaria may be significant [30]. The ethnopharmacological activities of this plant, including anti-inflammatory, antidiarrheal, antimicrobial, antioxidant, thrombolytic, hepatoprotective, analgesic, larvicidal, and insecticidal activities, have been reported in previous studies [25,26,31]. Importantly, the plant exhibits significant antioxidant and thrombolytic properties, which may be useful against these components of the SARS-CoV-2 pathogenesis. This study computationally investigated the roles played by compounds identified in L. zeylanica to combat SARS-CoV-2.

Collection and Identification of Plant Material
The aerial parts of L. zeylanica were collected from the forest area of Chittagong Hill Tract in November 2015, which was acknowledged by a prominent botanist from Bangladesh Forest Research Institute, Chittagong, Bangladesh. The Bangladesh National Herbarium has stored this plant for future reference with a voucher specimen for identification (accession no. BFRI-107).

Preparation of Plant Extract and Decoction Preparation
The aerial parts of the L. zeylanica plant were cleaned and cut into tiny pieces (0.4-0.5 mm), dried in air and ground (Moulinex Blender AK-241, Moulinex, France) into a powder (40-80 mesh, 355 g). The powder was immersed in 2 L of methanol for approximately 5 min, and the decoction was allowed to stand for 30 min before being filtered through Whatman No. 1 filter paper. At a temperature of less than 50 • C, filtrate from cheesecloth and Whatman filter paper No. 1 was condensed using a rotating evaporator (RE 200, Bibby Sterling Ltd., Staffordshire, UK). The extracts were placed into glass Petri dishes (90 × 15 mm, Pyrex, Germany) and allowed to air dry to evaporate the solvent completely, resulting in a final yield of 4.4%-5.6% w/w.

Gas Chromatography-Mass Spectrometry (GC-MS) Analysis
The aerial parts of L. zeylanica (methanol extract) were inspected with a mass spectrometer (TQ 8040, Shimadzu Corporation, Kyoto, Japan), using an electron-impact ionization process, combined with a gas chromatograph (GC-17A, Shimadzu Corporation), using a silica capillary (Rxi-5 ms; 0.25 m film, 30 m long and internal diameter 0.32 mm). The oven temperature was set at 70 • C (0 min); 10 • C, 150 • C (5 min); 12 • C, 200 • C (15 min); 12 • C, 220 • C (5 min), with a hold time of 10 min. The inlet temperature was 260 • C. A rate of 0.6 mL/min was used for the flux at a constant pressure of 90 kPa using helium gas. The temperature interface between the GC and MS was maintained at a constant 280 • C. The MS was performed in scan mode, with a range from 40 to 350 amu. The sample was injected at 1 µL, and the entire GC-MS process lasted for 50 min [32]. The results were compared against the National Institute of Standards and Technology (NIST) GC-MS library version 08-S for the identification of compounds in the peak areas.

Ligand Preparation
From the GC-MS analysis, 34 compounds (a total of 35 with the standard) were downloaded in .sdf two-dimensional (2D) format from the PubChem database (https: /pubchem.ncbi.nlm.nih.gov/, accessed on 12 March 2021). For the preparation of the ligand, the LigPrep tool (Maestro v 11.1) was used. The pH 7.0 ± 2.0 was used for the generation of ionization states of the compounds with Epik 2.2 (Force field: OPLS3) in Schrödinger ver.11.1. Up to 32 possible stereoisomers per ligand were retained.

Receptor Grid Generation and Glide Molecular Docking
The grid generation (Schrödinger ver.11.1) for the selected receptor was performed using the default parameters (Force field: OPLS3). Receptor grids were calculated for the prepared proteins for the observation of poses by various ligands bound within the active predicted site during the docking procedure. The van der Waals radius scaling factor and the partial atomic charge were 1.00 and 0.25, respectively. A cubic box of specific dimensions centered on the centroid of the active site residues was obtained for the receptor. The bounding box was set to 14 × 14 × 14 Å for docking experiments. Ligand docking was followed by the flexible standard precision (Schrödinger ver.11.1), and the docking score and the interactions of the ligand docking were recorded [35]. The results were represented as negative scores in kcal/mol, the final scoring was done on energy-minimized poses and shown as a Glide score. Discovery Studio (DS) version 4.5 was utilized to generate the 2D and 3D representations of the compounds [36]. The figures were generated using Microsoft PowerPoint 2019.

ADME/T Properties Analysis
The pharmacokinetic characteristics of all identified phytocompounds were evaluated and screened for drug candidacy using Lipinski's rule of five (RO5) [37]. According to Lipinski's RO5, a compound may exhibit optimal drug-like behavior if the selected parameters are within the specified limit and do not violate more than one of the following five criteria: molecular weight <500 g/mol; ≤5 hydrogen bond donors; ≤10 hydrogen bond acceptors; lipophilicity <5; and molar refractivity between 40 and 130. The web-based tool SwissADME, which is considered to be a convenient drug discovery tool, was used to analyze the druglikeness criteria of the detected biological compounds [38]. Compounds that pass Lipinski's RO5 can be considered suitable candidates for the development of new drug entities.

Molecular Dynamics Simulation
The dynamic motion of the drug-protein complex was evaluated through a molecular dynamics simulation study. The simulation study was conducted using the YASARA software (version 20.1.1) [39] with the aid of the AMBER14 force field [40]. The N3 inhibitor complex was used as a control system to be compared against the other three complexes. The periodic boundary condition was maintained, and a cubic simulation cell was created that was 20 Å larger than the biological systems in all cases. The NPT ensemble method was used and a Berendsen thermostat was applied to control the temperature of four systems. For the calculation of long-range electrostatic interactions, the particle mesh Ewald method [41] was applied, and the short-range van der Waals and Coulomb interactions were analyzed using a cutoff radius of 8 Å. The TIP3P or transferable intramolecular potential 3 points were used to add Na and Cl ions [42]. The total physiological conditions of the system were set to a temperature of 298 K, pH 7.4, and 0.9% NaCl. For the initial energy minimization process, the steepest descent gradient approach was used with a simulated annealing method. The molecular dynamics simulation trajectories were saved after every 100 ps using a time step of 1.25 fs [43]. The simulation study was conducted for 100 ns to analyze the root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), radius of gyration (Rg), solvent-accessible surface area (SASA), secondary structure, and the number of hydrogen bonds [44][45][46]. The molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method was applied to calculate the binding free energy. The YASARA macro file was edited for this calculation [47,48]. The 1000 trajectory files were considered for MM-PBSA calculation.

GC-MS Analysis
A total of 34 compounds were identified from the aerial parts of L. zeylanica using GC-MS, which are listed in Figure 1 and Table 1, along with their chemical compositions. The total ionic chromatogram (TIC) is shown in Figure 1. Thirty-four compounds were selected for molecular docking analyses because the specific biological activities of interest for these compounds have not yet been established.

Molecular Docking Study
The interactions between various identified compounds from L. zeylanica and the SARS-CoV-2 receptor (PDB ID: 6LU7) are presented in Table 2. Of the 34 compounds analyzed, 26 compounds interacted with the SARS-CoV-2 receptor, and the three compounds with the highest docking scores were identified as 11-oxa-dispiro  Figure S1. Importantly, the docking experiment delineated that standard inhibitor N3 showed the highest docking score (−7.013 kcal/mol) for the SARS-CoV-2 M pro compared to the other studied compounds.

GC-MS Analysis
A total of 34 compounds were identified from the aerial parts of L. zeylanica using GC-MS, which are listed in Figure 1 and Table 1, along with their chemical compositions. The total ionic chromatogram (TIC) is shown in Figure 1. Thirty-four compounds were selected for molecular docking analyses because the specific biological activities of interest for these compounds have not yet been established.

Molecular Dynamics Simulation
In this simulation study, lorazepam, 11-oxa-dispiro[4.0.4.1]undecan-1-ol, azetidin-2one-3, 3-dimethyl-4-(1-aminoethyl), and inhibitor N3 were denoted as D1, D2, D3, and control, respectively. The RMSD from Figure 6A illustrated that the control drug-protein complex had a higher RMSD trend compared with those of the other three complexes. Initially, none of the four complexes displayed large fluctuations, and they generally remained in a steady state. However, after 40 ns, all complexes had slightly higher and lower RMSD trends, which indicated complex flexibility. Eventually, the complexes returned to a steady trend again for the remainder of the simulation time, exhibiting rigid profiles. The degree of mobility in a biological system can be indicated by the Rg profile. Figure 6B shows that the control drug-protein complex had a lower Rg profile than the experimental drugs, indicating the compacted nature of the protein complex, whereas higher Rg values, which correlate with the repeated folding and unfolding protein behavior, were observed for the protein complexes containing both D1 and D2. By contrast, the complex containing D3 demonstrated initial stability until 40 ns, after which the Rg value increased, which may represent the loose packaging of the system. However, after 60 ns the D3 complex regained a steady Rg pattern, similar to those observed for the other complexes.
The surface area of the biological systems and their corresponding binding patterns with ligand molecules can be assessed through SASA analysis. The D3 complex showed an increasing SASA value until 40 ns, after which the SASA stabilized ( Figure 6C). This increasing trend in SASA represents protein expansion and comparatively loose binding. By contrast, D1, D2, and the control complex presented with stable SASA profiles and did not deviate. Therefore, these complexes experienced no changes in the surface area and formed more rigid profiles compared with the D3 complex. The quantitative measurement of hydrogen bonds in the drug-protein complex represents the constant nature and molecular recognition of the complexes. The D1 and D3 complexes formed more hydrogen bonds than the control and D2 complexes. More hydrogen bonds indicate an increasingly stable nature for the complex ( Figure 6D).

Discussion
With an increasing number of cases worldwide, the COVID-19 situation con worsen on a daily basis, especially in developed countries such as the USA. Th giousness of SARS-CoV-2 is above and beyond that of SARS-CoV or MERS-C In molecular modeling and computational drug design schemes, the PBSA system is a widely used solvation model for the estimation of the binding free energy of the drugprotein complex. Better binding and more compact interactions are indicated by higher binding energy values in the MM-PBSA calculations [44][45][46]. The reference control protein structure of the protein had more binding energy, which indicated a positive interaction pattern. The complex molecules D2 and D3 had similar binding energy patterns to that observed for the control complex, whereas D1 had slightly less binding energy than the control complex, as demonstrated in Figure 6F.

Discussion
With an increasing number of cases worldwide, the COVID-19 situation continues to worsen on a daily basis, especially in developed countries such as the USA. The contagiousness of SARS-CoV-2 is above and beyond that of SARS-CoV or MERS-CoV, two other members of the beta-coronavirus family. However, no specific treatment has been developed to treat this infectious disease thus far. As a result, the mortality rate continues to increase rapidly. Due to the absence of specific therapeutic drugs, current treatment approaches primarily involve symptom relief and supportive care. Several clinical trials have been conducted for several drug candidates, including remdesivir, hydroxychloroquine, and lopinavir/ritonavir [49]. However, these drugs have not yet amassed sufficient evidence to support their clinical applications. Scientists and researchers worldwide are working together to identify treatment strategies for this deadly coronavirus. Recently, a research study suggested a role for immunopathological considerations in the treatment of SARS-CoV-2 [50]. Wang et al. showed that human monoclonal antibodies could have a neutralizing effect against SARS-CoV-2, primarily by targeting the spike glycoprotein of the virus [51]. However, traditional drug development processes are tedious, and the development of suitable drug candidates can take as long as 15 years, which is not a feasible approach for identifying appropriate cures for COVID-19. Computer-aided drug design may represent a potential method for identifying lead compounds to treat SARS-CoV-2, which can result in both rapid and accurate results. A plethora of studies that have applied computational biology techniques have successfully predicted novel lead compounds for combating SARS-CoV-2 in addition to designing in silico epitope-based vaccines against SARS-CoV-2 [2,5,30,52,53]. However, these types of studies continue to require further wet-lab verification to be developed into effective therapeutic strategies for COVID-19.
Although drug research has indicated a high level of interest in the investigation of natural products, realistically, the synthesis and purification of vast arrays of new compounds represent a research bottleneck [74]. To compensate for the time-consuming and costly nature of new product development, high-throughput screening techniques have been developed as effective methods for the identification of new hit compounds [75]. However, these hit compounds often result in pharmacokinetic failures that result in elimination following Phase II clinical trials [76]. Therefore, current research works have begun to involve key investigations of pharmacokinetics and pharmacodynamics parameters, which are often referred to as drug-likeness properties [77]. Several rules have been established to facilitate the acceptance of promising molecules, and Lipinski's RO5 is the most commonly applied set of drug-likeness attributes [78]. In the current study, we analyzed the druglikeness properties of the identified compounds based on Lipinski's RO5 by utilizing the SWISS-ADME server. Our analysis indicated that most of the identified phytocompounds derived from the aerial parts of L. zeylanica followed Lipinski's RO5. Furthermore, it was reported that, the antiviral drugs approved by the FDA in the last five years involve a molecular weight of 513.97, hydrogen bond donors of 2.95, and hydrogen bond acceptors of 9.13 [23]. The targeted natural compounds are also within the reported criteria, indicating that these compounds may be evaluated as potential drug molecules.

Conclusions
We used a methodology for computer-aided drug design to detect potent SARS-CoV-2 M pro inhibitors following the L. zeylanica phytochemical analyses. Our research showed that several of the identified phytocompounds from L. zeylanica followed Lipinski's RO5 and were therefore assessed as promising therapeutic molecules. The catalytic residues M pro , Ser 144, Cys 145, and His41 were shown to be linked in postmolecular dynamic structures to the investigated phytochemicals. The molecular dynamics simulations conducted for the docked complexes revealed more insight into the rigidity and binding stability of these protein-drug complexes. Following additional investigations in biological laboratories, these compounds may possibly be developed into effective SARS-CoV-2 therapeutic candidates.