A Computational Study of Carbazole Alkaloids from Murraya koenigii as Potential SARS-CoV-2 Main Protease Inhibitors.

Despite COVID-19 vaccination, immune escape of new SARS-CoV-2 variants has created an urgent priority to identify additional antiviral drugs. Targeting main protease (Mpro) expressed by SARS-CoV-2 is a therapeutic strategy for drug development due to its prominent role in viral replication cycle. Leaves of Murraya koenigii are used in various traditional medicinal applications and this plant is known as a rich source of carbazole alkaloids. Thus, this computational study was designed to investigate the inhibitory potential of carbazole alkaloids from Murraya koenigii against Mpro. Molecular docking was initially used to determine the binding affinity and molecular interactions of carbazole alkaloids and the reference inhibitor (3WL) in the active site of SARS-CoV-2 Mpro (PDB ID: 6M2N).The top scoring compounds were further assessed for protein structure flexibility, physicochemical properties and drug-likeness, pharmacokinetic and toxicity (ADME/T) properties, antiviral activity, and pharmacophore modeling. Five carbazole alkaloids (koenigicine, mukonicine, o-methylmurrayamine A, koenine, and girinimbine) displayed a unique binding mechanism that shielded the catalytic dyad of Mpro with stronger binding affinities and molecular interactions than 3WL. Furthermore, the compounds with high affinity displayed favorable physicochemical and ADME/T properties that satisfied the criteria for oral bioavailability and druggability. The pharmacophore modeling study shows shared pharmacophoric features of those compounds for their biological interaction with Mpro. During the molecular dynamics simulation, the top docking complexes demonstrated precise stability except koenigicine. Therefore, mukonicine, o-methylmurrayamine A, koenine, and girinimbine may have the potential to restrict SARS-CoV-2 replication by inactivating the Mpro catalytic activity.

Introduction

COVID-19 has become the severest threat to the global public health in the twenty-first century causing more than 6 million mortalities worldwide [1]. The devastating pandemic has created unprecedented health, social, and economic crises all over the world. In the short span of 2 years, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has evolved to raise five variants of concern out of which Omicron has currently become the dominant variant all over the world [2]. A reduced efficiency of vaccines in antibody neutralization against Omicron variant has also been observed [3].

Currently, different vaccine programs are being used to combat with the pandemic. These vaccines are mainly four categories, viz. whole virus vaccines, viral vector vaccines, virus protein subunit vaccines, and mRNA vaccines that can trigger the immune system to produce antibodies against the virus [4]. However, COVID-19 vaccination is challenged by several factors such as vaccine hesitancy, failure of lifelong immunity, unequal distribution among countries, paucity of supply, and low efficacy for new variants [5]. Therefore, it is a pressing need to discover additional drug candidates for the treatment of COVID-19 infections.

Among the key targets of SARS-CoV-2 for antiviral therapy, main protease (Mpro) has attracted the interest of drug research efforts. The SARS-CoV-2 Mpro (SC2-Mpro) autoproteolytically cleaves the overlapping polyproteins (pp1a and pp1ab) to produce mature non-structural proteins (11 proteins) required for viral replication and transcription [6]. Intriguingly, the absence of a similar human protease to Mpro, as well as its unique cleavage site selectivity, has made it one of the primary therapeutic targets for anti SARS-CoV-2 drug discovery [6]. In comparison to the wild-type strain, the novel Omicron variant has one mutation in its major protease (P132H) (Supplementary Fig. S1). Nonetheless, it is not near the active site and may not interfere with the effectiveness of antiviral medications that target the active site. Paxlovid marketed by Pfizer Inc. is the only drug licensed for the use as an oral antiviral drug against Mpro protein under emergency use authorization (EUA) [7].

The absence of effective therapeutic drugs has boosted the desire for natural alternatives to treat COVID-19. A substantial amount of research has reported the antiviral activity of plant derived compounds against RNA viruses [8]. Murraya koenigii, often known as the “Curry leaf tree,” is a tropical and subtropical plant that is extensively grown across the world. In Sri Lanka, its leaves are widely used for cooking and preparing leafy porridge along with other medicinal applications [9].

The plant has been used to treat various diseases in Ayurveda system utilizing its wide spectrum of pharmacological activities such as antioxidant, antidiabetic, anticancer, anti-inflammatory, hepatoprotective, nephroprotective, cardioprotective, neuroprotective, and antimicrobial activities [10]. However, the antiviral activity of Murraya koenigii has been less explored. Previous research reports have shown that the leaves, roots, and bark of this plant are high in carbazole alkaloids [10]. These alkaloids have a variety of biological actions, including antioxidant, anti-diabetic, anti-inflammatory, anticancer, and neuroprotective properties [10]. The carbazole alkaloids include a carbazole scaffold, which is a favored scaffold that is significant in drug development. Many drug compounds containing a carbazolic core have been discovered, and some have demonstrated antiviral action [11]. The goal of this computational study was to identify potent SC2-Mpro inhibitors among Murraya koenigii carbazole alkaloids (8 compounds), which may be useful in developing antiviral drugs for COVID-19.

Materials and Methods

Preparation of Structures for Molecular Docking

The crystal structures of viral protein, SC2-Mpro (PDB ID: 6M2N, 2.20 A° resolution, PDB ID: 7U29, 2.09 A° resolution; PDB ID: 7TLL, 1.63 A° resolution), were procured from the Protein Data Bank database (https://www.rcsb.org/). Except one crystallographic water molecule (532H2O), all the rest of the water molecules, ligands, and heteroatoms were removed from the crystal structure of PDB ID: 6M2N. All water molecules, ligands, and heteroatoms were removed from the remaining crystal structures (PDB ID: 7U29 and 7TLL). All structures were adjusted to the physiological protonation state using UCSF Chimera software [12]. Single chain A of all SC2-Mpro structures was selected and saved in PDB format for further molecular simulations.

Previously reported eight natural carbazole alkaloids from Murraya koenigii were compiled from the literature [10]. They were imported into Avogadro v.1.2.0 software [13] using PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Then, all carbazole alkaloids were adjusted to physiological protonation states followed by energy minimization in the universal force field (UFF) with conjugate gradient algorithm for 500 steps using Avogadro v.1.2.0 software. The co-crystal ligand, 3WL (baicalein) was used as the reference inhibitor which was taken from the respective PDB entry of SC2-Mpro (PDB ID: 6M2N). It was adjusted to physiological protonation state by Avogadro v.1.2.0 software.

Identification of Active Site Residues

The corresponding publication of the SC2-Mpro protein structure (PDB ID: 6M2N) [14] was referred to identify the crucial active site residues. This was crosschecked by analyzing the 5 A° zone of the co-crystal ligand, 3WL bound at the active site of the SC2-Mpro using UCSF Chimera software.

Molecular Docking and Molecular Interaction Studies

Molecular docking was executed to identify the binding affinity and interactions between wild-type SC2-Mpro (PDB ID: 6M2N) and each compound using AutoDock 4.2 software [15]. Re-docking of the co-crystal ligand, 3WL was exploited to validate the docking experiment. It was confirmed by computing the root-mean square deviation (RMSD) after superimposing native crystal pose with the re-dock pose. The grid box was placed enclosing the active site residues of the SC2-Mpro with the dimensions of 40 A° whereas center grid box parameters were positioned at − 32.981, − 65.109, and 41.712 in X, Y, and Z axes, respectively. During the docking process, the Kollman and Gastieger charges were added to the protein and ligands, respectively. Molecular docking calculations were performed in Lamarckian genetic algorithm (LGA). Rigid protein structure and flexible ligands were used throughout the molecular docking process. Ligand efficiency was calculated by dividing the binding affinity (-∆G) by the number of non-hydrogen atoms in the ligand (n-NHA). The binding energy of the co-crystal ligand (3WL) was set as the cutoff value to screen the top scoring phytochemicals for further analysis.

Moreover, the findings of AutoDock 4.2 were validated applying AutoDock Vina software [41]. A grid box of 15 × 15 × 15 A° dimensions was prepared to target the active site residues and spacing was kept at 1 A°. The maximum energy difference between the best and worst binding modes was set at 2 kcal/mol with an exhaustiveness value of 100.

BIOVIA Discovery Studio Visualizer Client 2020 [16] was used to conduct a detailed study of molecular interactions between the wild-type SC2-Mpro and compounds. The two dimensional (2D) and three dimensional (3D) interaction plots of the top docking complexes identified from AutoDock 4.2 were visualized using BIOVIA Discovery Studio Visualizer Client 2020.

Further top scoring phytochemicals were docked with the SC2-Mpro of the previously evolved variants (Alpha, Beta, and Gamma) and the current variant (Omicron) using PDB ID: 7U29 and 7TLL structures respectively. Alpha, Beta, and Gamma (K90R) and Omicron (P132H) variants carry mutations in the SC2-Mpro. However, the Delta variant does not carry any mutation in the Mpro sequence and carries the same sequence as the wild-type.

Molecular docking was carried out between the SC2-Mpro of variants and the top scoring phytochemicals using AutoDock Vina software. With regard to PDB ID: 7U29, the center grid box parameters were positioned at − 5.78, 2.524, and 10.702 in X, Y, and Z axes, while the center grid box parameters were positioned at − 3.107, 0.755, and 13.946 in X, Y, and Z axes for PDB ID: 7TLL, respectively. All the other parameters were set as the previously mentioned under the procedure for molecular docking analysis using AutoDock Vina.

Prediction of In Silico Physico-chemical Properties and Drug-Likeness

Molinspiration web tool (https://www.molinspiration.com/) [18] was employed to compute several physico-chemical properties in silico and to predict drug-likeness for the top scoring phytochemicals. The drug-likeness was predicted using Lipinski’s rule [19] and Veber’s rule [20].

Molecular Dynamics Simulations

The best-docked protein–ligand complexes based on AutoDock 4.2 were subjected to molecular dynamics (MD) simulation for 100 ns utilizing GROMACS simulation package based WebGro server (https://simlab.uams.edu/). The simulation of protein–ligand complexes were performed using GROMOS96 43a1 force field. The ligand topology files were generated by PRODRG web tool (http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg) [17]. This system was solvated using SPC water model in a triclinic box and neutralized adding 0.15 M counter ions (Na + /Cl −). Further steepest descent algorithm (5000 steps) was applied to minimize energy of the system. Finally, MD simulation was run using Leap-frog integrator with NVT/NPT equilibration under temperature (300 K) and pressure (1.0 bar). Approximate number of frame per simulation was set to 1000. The MD trajectories were used to analyze the stability and conformational changes of the protein–ligand complex using different parameters, namely root mean square deviation (RMSD), root mean square fluctuation (RMSF), and radius of gyration (Rg).

In Silico ADME/T Predictions

Pharmacokinetic properties of the top hit phytochemicals, related to absorption, distribution, metabolism, and excretion (ADME), were predicted through the pkCSM web server (http://biosig.unimelb.edu.au/pkcsm/) [21], while toxicity (T) analysis was performed using ProTox-II web server (https://tox-new.charite.de/) [22].

In Silico Prediction of Antiviral Activity

Prediction of antiviral activity for the selected phytochemicals was ascertained by AVCpred web server (http://crdd.osdd.net/servers/avcpred/) [23]. The web server was designed based on quantitative structure–activity relationship (QSAR) models that were built using known antiviral compounds against human immunodeficiency virus (HIV), hepatitis C virus (HCV), hepatitis B virus (HBV), human herpesvirus (HHV), and 26 other important viruses with experimentally validated percentage inhibition from ChEMBL database.

Pharmacophoric Feature Analysis and Pharmacophore Model Building

The mol2 format of top 5 phytochemicals was subjected to pharmacophoric feature analysis and pharmacophore model building through PharmaGist web server (https://bioinfo3d.cs.tau.ac.il/PharmaGist/) [24]. The merged pharmacophore model was built based on the top score after structural alignment of top 5 phytochemicals. The phytochemical with the lowest binding energy was selected as the key molecule and all other top scoring phytochemicals were flexibly aligned onto it in order to build the merged pharmacophore model. The output files obtained from PharmaGist web server were visualized using ZINCPharmer web server (http://zincpharmer.csb.pitt.edu/) [25].

Results

The SC2-Mpro protein consists of three domains: domain I (containing residues 8–101), domain II (containing residues 102–184), and domain III (containing residues 201–303). The active site residues lie on the surface of the protein between the domains I and II. They are mainly composed of a catalytic dyad (His41 and Cys145) and several binding site residues (Leu27, Cys44, Asp48, Met49, Pro52, Tyr54, Leu141, Asn142, Gly143, Ser144, His163, His164, Met165, Glu166, Leu167, Asp187, Arg188, and Gln189). The aforementioned active site residues were used to design the grid box for further molecular docking studies.

Natural carbazole alkaloids along with 3WL were docked into the active site of the target protein using molecular docking experiment in order to identify the best docking pose with adequate molecular interactions.

The re-docking experimental results of AutoDock 4.2 and AutoDock Vina produced a root mean square deviation (RMSD) value of 0.219 A° and 0.268 A°, respectively, after superimposition of co-crystal ligand (3WL) pose with its re-docking ligand pose (Fig. 1 a-i and a-ii). This confirmed the ability of the docking experiment to reproduce native poses.

Fig. 1

a Superimposition of native ligand pose (yellow) and re-docking ligand pose (green) of (i) AutoDock 4.2 and (ii) re-docking ligand pose (cornflower blue) of AutoDock Vina. b Overlay of top scoring compounds and 3WL in the binding pocket 3WL, yellow; koenigicine, blue; mukonicine, purple; o-methylmurrayamine A, green; koenine, orange; girinimbine, gray

It was discovered that all the compounds interact with the active site of the protein. According to the AutoDock 4.2 results, the binding energies of carbazole alkaloids ranged from − 7.26 to − 6.57 kcal/mol, whereas the inhibition constants varied between 4.78 and 15.25 μM. The co-crystal ligand (3WL) which was also considered the reference inhibitor showed a binding energy of − 6.84 kcal/mol with an inhibition constant (Ki) of 9.70 μM. All the carbazole alkaloids exhibited high ligand efficiencies ≥ 0.30 except koenigine which showed a value of 0.29. Among the selected carbazole alkaloids, five compounds (koenigicine, mukonicine, o-methylmurrayamine A, koenine, and girinimbine) recorded lower values for binding energy (higher docking score) than the reference inhibitor (3WL), while remaining three compounds (murrayacine, koenigine, and koenimbine) showed higher values. The computed values for the docking parameters including binding energy, inhibition constant, and ligand efficiency of all compounds are listed in Table 1. On the other hand, AutoDock Vina resulted binding energies of carbazole alkaloids ranged from − 7.8 to − 6.8 kcal/mol, whereas the co-crystal ligand, 3WL (reference inhibitor), showed a binding energy of − 7.2 kcal/mol (Table 1). The molecular docking results of SC2-Mpro in different variants are tabulated in Table 2.

Table 1

Docking results and 2D structures of the compounds of Murraya koenigii used for this study

*Ligand Binding Energy, **Estimated Inhibition Constant, ***Ligand Efficiency

N/A – Not Applicable

Table 2

AutoDock Vina results of the compounds of Murraya koenigii against SARS-CoV-2 Mpro variants (Alpha, Beta, Gamma, and Omicron)

		Ligand binding Energy (kcal/mol)
No	Ligand Name	Alpha, Beta, Gamma (K90R)	Omicron (P132H)
		(PDB ID: 7U29)	(PDB ID: 7TLL)
1	Koenigicine	− 6.6	− 6.9
2	Mukonicine	− 7.2	− 7.5
3	O-Methylmurrayamine A	− 6.8	− 7.1
4	Koenine	− 7.3	− 7.4
5	Girinimbine	− 7.5	− 7.5
6	3WL (reference)	− 7.5	− 7.5

Moreover, BIOVIA Discovery Studio Visualizer was employed to study the molecular interactions of the best dock poses in details. Inside the active site, all nine compounds were stabilized by various non-covalent interactions such as conventional hydrogen, carbon hydrogen, pi-sigma, alkyl, pi-alkyl, pi-pi stacked, and pi-sulfur contacts. These non-covalent interactions were mediated through 19 residues, namely Leu27, His41, Cys44, Asp48, Met49, Pro52, Tyr54, Leu141, Asn142, Gly143, Ser144, His163, His164, Met165, Glu166, Leu167, Asp187, Arg188, and Gln189. The recurrent residues of conventional hydrogen bond formation are Gly143 (5/9), Glu166 (3/9), His41 (2/9), and Asn142 (2/9). The residues, Gln189, Cys145, and Ser144, formed only one conventional hydrogen bond. The detailed analysis of conventional hydrogen bond interactions of the best dock poses is shown in Table 3. It was noticeable that His41 was able to form at least one of the non-covalent interactions with all compounds. The heat maps of protein–ligand molecular interaction studies based on AutoDock 4.2 and AutoDock Vina are illustrated in Fig. 2 a and b separately.

Table 3

Conventional hydrogen bond interaction analysis of the docking complexes

No	Docking complex	nCHBa	CHBIb		Distance (A°)
			Protein	Ligand
1	Mpro—Koenigicine	1	A: Gly143: HN	UNL1: O	2.51
2	Mpro—Mukonicine	2	A: Gly143: HN	UNL1: O	2.69
			A: His41: ND1	UNL1: H	2.92
3	Mpro—O-Methylmurray	2	A: Asn142: HD22	UNL1: O	2.85
	-amine A		A: Gln189: HN	UNL1: O	3.07
4	Mpro – Koenine	2	A: Glu166: HN	UNL1: O	2.53
			A: Asp187: O	UNL1: H	1.82
5	Mpro – Girinimbine	2	A: Asn142: HD22	UNL1: O	3.04
			A: His41: O	UNL1: H	3.08
6	Mpro – 3WL*	2	A: Gly143: HN	A: 3WL401: O3	2.33
			A: Glu166: HN	A: 3WL401: O	2.04
7	Mpro – Murrayacine	4	A: Gly143: HN	UNL1: O	2.34
			A: Gly143: HN	UNL1: O	2.19
			A: Ser144: HN	UNL1: O	2.73
			A: Cys145: HN	UNL1: O	2.37
8	Mpro – Koenigine	2	A: Gly143: HN	UNL1: O	2.66
			A: Asp187: O	UNL1: H	2.21
9	Mpro – Koenimbine	1	A: Glu166: HN	UNL1: O	2.72

aNumber of conventional hydrogen bonds (nCHB)

bConventional hydrogen bond interaction (CHBI); donor atoms, bold; acceptor atoms, italics

*3WL forms one water hydrogen bond with a crystallographic water molecule (HOH532)

The distance of the bond is 2.18 A°

Fig. 2

Protein–ligand molecular interaction studies of all compounds using a AutoDock 4.2 and b AutoDock Vina (two catalytic residues are shown in red color and numerical values inside the diagram denote the number of interactions)

The oxygen atom (heteroatom) of 2H-pyran ring of koenigicine formed a hydrogen bond with Gly143. The amine group (-NH) in the carbazole scaffold of mukonicine formed one hydrogen bond contact with the catalytic site residue, His41, and another hydrogen bond with Gly143 through the heteroatom of 2H-pyran ring. The heteroatom of 2H-pyran ring and methoxy group (-O-CH3) of carbazole scaffold of o-methylmurrayamine A formed one hydrogen bond each with Asn142 and Gln189, respectively. Koenine formed two hydrogen bonds with Glu166 and Asp187 of SC2-Mpro interacting with the heteroatom of the 2H-pyran ring and hydroxyl group of carbazole scaffold, respectively. The heteroatom of 2H-pyran ring of girinimbine formed a hydrogen bond with Asn142, while the amine group (-NH) in the carbazole scaffold interacted with His41 to form another hydrogen bond. The heteroatom of 2H-pyran ring of murrayacine formed a hydrogen bond with Gly143 and the carbonyl oxygen of aldehyde group established three hydrogen bonds with Cys145, Ser144, and Gly143. Koenigine formed two hydrogen bond contacts with Gly143 and Asp187 through the heteroatom of 2H-pyran ring and hydroxyl group, respectively. The heteroatom of 2H-pyran ring of koenimbine formed only one hydrogen bond with Glu166. The reference inhibitor established two hydrogen bonds with Glu166 and Gly143 through the carbonyl oxygen and hydroxyl group (6-OH), while another hydroxyl group (5-OH) formed one water hydrogen bond with a crystallographic water molecule (532H2O). The 2D and 3D molecular interaction diagrams of the best docking poses of nine compounds are depicted in Fig. 3.

Fig. 3

The 3D and 2D molecular interaction diagrams of the best docking poses of compounds

Among all carbazole alkaloids, koenigicine, mukonicine, o-methylmurrayamine A, koenine, and girinimbine which exhibited better binding affinities as well as stronger molecular interactions than 3WL were selected for further molecular dynamics simulations along with 3WL.

Numerous physico-chemical descriptors of top five phytochemicals and 3WL were assessed for druggability using Molinspiration webserver. Further drug-likeness rules such as Lipinski’s rule of five (LRo5), molecular weight (MW) =  < 500, number of hydrogen bond acceptors (nHA) =  < 10, log of the octanol/ water partition coefficient (LogP) =  < 5, number of hydrogen bond donors (nHD) =  < 5 and Veber’s rule; number of rotatable bonds (nRot) =  < 10, and topological polar surface area (TPSA) =  < 140 were applied to identify orally active compounds. The predicted descriptors of top five phytochemicals such as MW, volume, nHA, nHD, nRot, nNHA, TPSA, and LogP lied within the ranges of 263.34–323.39 g mol−1 (MW), 248.76–299.85 A°3 (volume), 2–4 (nHA), 1–2 (nHD), 0–2 (nRot), 20–24 (nNHA), 25.02–45.25 A° (TPSA), and 4.52–5.05 (LogP), respectively (Table 4). Both koenigicine and mukonicine displayed comparable results for physicochemical descriptors except LogP values which were slightly changed. Three compounds (mukonicine, o-methylmurrayamine A, and girinimbine) surpassed the maximum limit of LogP value (5), while remaining compounds showed low values. As one violation of the descriptors is permitted for LRo5, all top hits satisfied the criteria for LRo5 and Veber’s rule. Hence, top five phytochemicals showed favorable drug-likeness properties.

Table 4

Physicochemical properties and drug-likeness properties of top hit compounds

No	Description	MWa	Volumeb	nHAc	nHDd	nRote	nHetf	TPSAg	LogPh	LRi	VRj
		(g/ mol)	(A°3)					(A°)
1	Koenigicine	323.150	339.163	4	1	2	4	43.480	5.209	Yes	Yes
2	Mukonicine	323.150	339.163	4	1	2	4	43.480	5.564	Yes	Yes
3	O-Methylmurraya	293.140	313.076	3	1	1	3	34.250	5.600	Yes	Yes
	mine A
4	Koenine	279.130	295.780	3	2	0	3	45.250	4.999	Yes	Yes
5	Girinimbine	263.130	286.990	2	1	0	2	25.020	5.565	Yes	Yes
6	3WL	270.050	265.186	5	3	1	5	90.900	3.215	Yes	Yes

aMolecular Weight

bVan der Waals Volume

cNumber of hydrogen bond acceptors

dNumber of hydrogen bond donors

eNumber of rotatable bonds

fNumber of heteroatoms

gTopological polar surface area

hLog of the octanol/ water partition coefficient

iLipinski’s rule of five (MW =  < 500, nHA =  < 10, LogP =  < 5, nHD =  < 5) and jVeber’s rule (nRot =  < 10, TPSA =  < 140)

Molecular dynamics simulations were exploited to analyze the behavior of protein’s stability after binding of phytochemicals to its active site.

The RMSD trajectory of the selected dock complexes is illustrated in Fig. 4a.The koenigicine complex showed high fluctuations after 35 ns till the end of simulation. All the other complexes attained a stable state after 40 ns.

Fig. 4

Molecular dynamics simulation plots of top scoring docking complexes, a RMSD plot, b RMSF plot, and c Rg plot

The RMSF values obtained for top scoring dock complexes and 3WL are plotted in Fig. 4b. Higher fluctuations were noticed in the C and N termini of all complexes. Less number of residues (79 of 306) of top five dock complexes and 3WL were observed in highly flexible region (> 3 A°). Out of 79 highly flexible residues, the majority belonged to the coil structures (Fig. 5). The percentages of highly flexible residues in each top dock complex were 16.67%( for koenigicine), 3.92% (for mukonicine), 13.40% (for o-methylmurrayamine A), 3.92% (for koenine), 3.27% (for girinimbine), and 10.78% (for 3WL) (Supplementary Table S1).

Fig. 5

Protein secondary structure analysis of SARS-CoV-2 M.pro (PDB ID: 6M2N) (the image was generated using PSIPRED web tool)

Initially, Rg values (Fig. 4c) of koenigicine complex showed significant variation compared to other five complexes till it reached equilibrium at around 40 ns. After that, Rg values of all six complexes were consistent around 2.1 nm till the end of the simulation time period. The calculated average values of Rg were found to be 2.13 nm, 2.13 nm, 2.13 nm, 2.12 nm, 2.14 nm, and 2.12 nm for 3WL complex, koenine complex, o-methylmurrayamine A complex, mukonicine complex, koenigicine complex, and girinimbine complex, respectively.

ADME/T assessments of the top five phytochemicals were determined using pkCSM and ProTox-II web servers. They are listed in Table 5. The top five phytochemicals were moderately soluble in the water and had almost 90% human intestinal absorption. Among the top five compounds, girinimbine exhibited the greatest value (0.798) for steady-state volume of distribution followed by o-methylmurrayamine A (0.624), koenigicine (0.6), mukonicine (0.516), and koenine (0.488). Mukonicine, o-methylmurrayamine A and koenine acted as a CYP3A4 substrate but inhibited CYP1A2 and 2C19. On the other hand, koenigicine and girinimbine inhibited all three CYP1A2, 2C9, and 2C19 enzymes but acted as a substrate for CYP3A4 enzyme. All the top five compounds showed moderate values for total clearance close to the minimal limit (0.301). No harmful toxicity was observed for koenine and girinimbine, while koenigicine, mukonicine, and o-methylmurrayamine A indicated carcinogenicity and mutagenicity.

Table 5

In silico pharmacokinetics and toxicity predictions of top hit compounds

No	Description	Absorption		Distribution	Metabolism		Excretion		Toxicity
		WSa	HIAb	VDssc	SUBd	INHe	TCf	HEPg	CARh	MUTi	CYTj
1	Koenigicine	-4.933	93.217	0.6	CYP3A4	CYP1A2	0.484	No	Yes	Yes	No
						2C19, 2C9
2	Mukonicine	-4.923	92.414	0.516	CYP3A4	CYP1A2, 2C19	0.514	No	Yes	Yes	No
3	O-Methylmurraya	-4.949	92.521	0.624	CYP3A4	CYP1A2, 2C19	0.483	No	Yes	Yes	No
4	Koenine	-4.516	89.944	0.488	CYP3A4	CYP1A2, 2C19	0.404	No	No	No	No
5	Girinimbine	-5.048	93.43	0.798	CYP3A4	CYP1A2, 2C19	0.381	No	No	No	No
						2C9
6	3WL	-3.302	94.268	-0.004		CYP1A2, 2C9	0.252	No	Yes	Yes	No

aWater solubility (WS)[log(mol/L)]

bHuman intestinal absorption (HIA)[%]

cSteady state volume of distribution (VDss) [log (L/kg)]

dSubstrate (SUB)

eInhibitor (INH)

fTotal clearance (TC) [log (ml/min/kg)]

gHepatotoxicity (HEP)

hCarcinogenicity(CAR)

iMutagenicity(MUT)

jCytotoxicity(CYT)

Accepted range of values: HIA %—> 30% / Water solubility [log (mol/L)]—insoluble < -10 < poorly soluble < -6 < moderately soluble < -

4 < soluble < -2 < very soluble < 0 < highly soluble / VDss [log (L/kg)]—high > 0.45, low < -0.15 / Total Clearance

[log (ml/min/kg)]—high > 1.176, low < 0.301

In Silico Antiviral Activity Predictions

Due to lack of experimental results for antiviral activity of lead compounds, AVCpred online tool was used to calculate in silico antiviral inhibition percentage of them. Koenine showed the highest percentage of antiviral inhibition of 57.944% on all tested general viruses. All top hits exhibited higher percentages of antiviral inhibition against HCV, HHV, and HIV compared to the 3WL. In contrast, 3WL showed higher percentage of antiviral inhibition against HBV than any of the top lead compounds. All the top hits expressed a considerable percentage of antiviral inhibition with all tested general viruses (Supplementary Table S2), HCV, HHV, and HIV compared to HBV. The predicted percentages of antiviral inhibition of each top hit compound and 3WL are tabulated in Table 6.

Table 6

Antiviral inhibition percentage prediction of top hit compounds

No	Ligand Name	Generala	HBVb	HCVc	HHVd	HIVe
1	Koenigicine	31.082	19.865	63.502	59.28	57.195
2	Mukonicine	33.428	18.748	59.145	53.232	55.049
3	O-Methylmurrayamine A	47.197	14.227	50.78	20.737	65.621
4	Koenine	57.944	13.245	57.802	49.038	55.525
5	Girinimbine	42.383	19.925	45.914	31.702	67.798
6	3WL (control)	48.864	31.751	37.526	20.068	33.709

aThe general dataset is comprised of 26 viruses including SARS coronavirus

bHepatitis B virus (HBV)

cHepatitis C virus (HCV)

dHuman herpesvirus (HHV)

eHuman immunodeficiency virus (HIV)

Analysis of Pharmacophoric Features and Pharmacophore Model Building

The 3D pharmacophoric features of the top scoring phytochemicals were elucidated using the PharmaGist webserver. The merged pharmacophore model was built considering koenigicine as the reference compound. The best pharmacophore model was selected based on the highest score and maximum features after flexible, structural alignment of top five phytochemicals. Several pharmacophoric features (pharmacophores), namely aromatic ring, hydrophobic nature, hydrogen bond donor, hydrogen bond acceptor, positive ionizable, and negative ionizable, were detected through the investigation. The maximum number of pharmacophoric features was found to be 14 for koenigicine and o-methylmurrayamine A, while the minimum number was recorded as 10 for girinimbine. The best pharmacophore model composed of 3 aromatic rings, 4 hydrophobic groups, 1 donor, 1 acceptor, and 1 positive ionizable group. The recognized pharmacophoric features of each lead compound and the best pharmacophore model are shown in Table 7 and Fig. 6.

Table 7

Three dimensional pharmacophoric features of the lead phytochemicals and the best pharmacophore model (overlapping top five hits)

No	Molecules	Atoms	Features	Spatial Features	Aromatic	Hydrophobic	Donors	Acceptors	N	P
1	A	45	14	13	3	6	1	3	0	1
2	B	38	12	10	3	4	2	2	0	1
3	C	45	14	13	3	6	1	3	0	1
4	D	41	12	11	3	5	1	2	0	1
5	E	37	10	9	3	4	1	1	0	1
6	A, B, C, D, E*	–-	10	9	3	4	1	1	0	1

A, koenigicine; B, mukonicine; C, o-methylmurrayamine A; D, koenine; E, girinimbine; N, negatives; P, positives

*Koenigicine molecule was used as the reference compound to build the topmost scoring pharmacophore model by overlapping the top five hits

Fig. 6

Three dimensional pharmacophoric features of the lead phytochemicals and the best pharmacophore model

Discussion

Computer-aided drug design techniques such as molecular docking have been extensively implemented by researchers for rapid recognition of potent inhibitors against viral targets of SARS-CoV-2 [26-28]. Out of several viral targets of SARS-CoV-2, Mpro plays a vital role in viral replication through proteolytic processing of the viral replicase polyproteins 1a and 1ab. Therefore, inhibition of the activity of Mpro is a therapeutic strategy to design potent inhibitors against COVID-19.

Presently, a new variant of SARS-CoV-2 named Omicron has dominated the world over Wuhan strain. Therefore, a pairwise sequence alignment of Mpro proteins’ sequences was conducted to identify any changes between the wild-type and Omicron variant (Supplementary Fig. 1). Interestingly, one mutation (P132H) was recognized in the Mpro of Omicron variant. However, the original active site of the protein is not affected by this mutation; thus, a potent inhibitor can interact with the active site of Mpro of Omicron variant in the same way as it does to wild-type Mpro without reducing its binding efficacy.

In the present study, SC2-Mpro (Wuhan strain) was used as a viral target to identify potent inhibitors from eight carbazole alkaloids present in Murraya koenigii through molecular docking simulation. Based on AutoDock 4.2 docking results, koenigicine (− 7.26 kcal/mol), mukonicine (− 7.14 kcal/mol), o-methylmurrayamine A (− 7.03 kcal/mol), koenine (− 6.92 kcal/mol), and girinimbine (− 6.89 kcal/mol) showed lower binding energies thus stronger binding affinity towards the active site of the SC2-Mpro than the reference inhibitor (3WL). It was found that 3WL showed the dock pose with binding energy of − 6.84 kcal/mol with SC2-Mpro closely followed by murrayacine, koenigine, and koenimbine with binding energies of − 6.82 kcal/mol, − 6.68 kcal/mol, and − 6.57 kcal/mol, respectively. All the selected carbazole alkaloids including 3WL exhibited ligand efficiency value ≥ 0.29 which is considered a good starting point for lead optimization [29].

Compared to 3WL, carbazole compounds exhibited significant interactions within the active site of the wild-type SC2-Mpro with 19 key residues (Leu27, His41, Cys44, Asp48, Met49, Pro52, Tyr54, Leu141, Asn142, Gly143, Ser144, His163, His164, Met165, Glu166, Leu167, Asp187, Arg188, and Gln189) out of which 8 residues (His41, Cys145, Asn142, Gly143, Ser144, Glu166, Asp187, and Gln189) participated in conventional hydrogen bond contacts.

The 2H-pyran ring of the selected carbazole alkaloids played a crucial role in forming conventional hydrogen bonds with the active site residues of SC2-Mpro. The oxygen atom of the 2H-pyran ring interacted to form one hydrogen bond each with Gly143 (koenigicine/mukonicine/murrayacine/koenigine), Asn142 (o-methylmurrayamine A/girinimbine), and Glu166 (koenine/koenimbine).

According to the bond interaction statistics, His41 interacted with all of the compounds by multiple non-covalent bond interactions, whereas Cys145 connected with all phytochemicals except koenimbine via multiple non-covalent bond interactions. Both of these residues are responsible for the catalytic activity of SC2-Mpro. As reported in a previous research [26], the formation of a conventional hydrogen bond contact or other non-covalent bond interactions with both Cys145 and His41 residues or any of these residues suggests the inhibitory potential of chosen compounds inside the active site of SC2-Mpro. Notably, koenine and koenimbine formed conventional hydrogen bond with Glu166 residue which is critical for the dimerization of SC2-Mpro monomers. This homodimer plays an essential role in the SC2-Mpro’s catalytic activity as the monomeric enzyme is not catalytically active [30]. As a result, any contact with Glu166 can obstruct the formation of a functional homodimer, which may eventually lead to a catalytically inactive SC2-Mpro protein [26]. Additionally, molecular interaction study unveiled that Gly143 in SC2-Mpro forms the highest number of hydrogen bonds (6) with the compounds (5/9) turning to an attractive site to form hydrogen bonds. This was in agreement with a previous study conducted using 137 crystal complexes of inhibitor bound SC2-Mpro [31].

Five carbazole compounds possessing lower binding energies and stronger molecular interactions than 3WL were identified for the further computational studies. All top five phytochemicals occupied in the same orientation in the vicinity of the catalytic dyad as the reference inhibitor.

Unlike other known SC2-Mpro inhibitors, the top scoring carbazole compounds showed a unique mode of interaction with SC2-Mpro. They acted as a shield in front of the catalytic site residues to block substrate reaching to the catalytic dyad (Fig. 1b-ii).This distinctive binding mode of inhibition along with lower binding energies and stronger molecular interactions makes the top five phytochemicals to be considered for further in silico analysis.

It was also observed that all top five phytochemicals including 3WL exhibited strong binding affinity towards the SC2-Mpro of different variants showing their efficacy to be developed as drugs against Alpha, Beta, Gamma, and Omicron variants. The wild-type and Delta variants have the same sequence for SC2-Mpro; thus, the top scoring phytochemicals can be expected to interact with the Mpro of Delta variant in the same way as it does to wild-type Mpro without any changes in the binding efficacy.

Furthermore, all top five phytochemicals were subjected to examine the physicochemical properties and drug-likeness. The physico-chemical descriptors of the top five phytochemicals were within the acceptable range satisfying both LR05 and Veber’s rules. This suggests their feasibility to be developed as an orally active antiviral drug.

Molecular dynamics simulations were conducted to understand the protein structure flexibility of each dock complex (top five) upon phytochemical binding.

The RMSD is often used to evaluate the stability of a ligand bound protein structure during MD simulation. The conformational stability of the protein–ligand complex will be higher if the RMSD value of complex is lower. It was seen through the RMSD trajectory that koenigicine complex exhibited higher fluctuations than the other five complexes hence showed the least stability. The RMSD values for the remaining five complexes (3WL, koenine, o-methylmurrayamine A, mukonicine, and girinimbine) fell within the acceptable range thus showed the formation of stable complexes with the SC2-Mpro.

The flexibility of individual amino acid was computed using RMSF calculations. Generally, larger RMSF values demonstrate high flexibility and low stability and vice versa. Higher peaks were observed at the C and N termini in the complexes indicating low stability due to high fluctuations. Totally, 79 residues showed high flexibility. They mainly appeared in the coil structures of the protein (Fig. 5). This reflects that alpha-helix and beta-sheet structures of the protein show mild fluctuations thus more rigid than the coils. Interestingly, none of the highly flexible residues associated in molecular interactions with mukonicine, koenine, and girinimbine. But four of highly flexible residues participated in molecular interactions with o-methylmurrayamine A (Asp48/Met49/Pro52/Gln189), while one highly flexible residue formed interactions with koenigicine (Gln189). Another highly flexible residue also involved in molecular interactions with 3WL. Koenigicine complex showed the maximum fluctuation of protein residues (16.67%) in comparison with the remaining five complexes (3WL, koenine, o-methylmurrayamine A, mukonicine, and girinimbine) thereby indicating the least stability among them. Therefore, all complexes of 3WL, koenine, o-methylmurrayamine A, mukonicine, and girinimbine demonstrated the precise rigidity of the protein structures upon attachment of compounds [32].

The Rg plot denotes the compactness of a ligand bound protein during simulation time. Lower protein compactness is indicated by a higher Rg value, and vice versa. Therefore, high Rg value indicates high flexibility in protein and low stability. Initially, koenigicine complex varied from low compactness to more compact till 40 ns and continued the constant compactness with relatively low Rg value (2.14 nm) for remaining 40 ns of simulation. However, all the other complexes exhibited high compactness (with low Rg values) indicating high stability throughout the simulation time.

Evaluation of pharmacokinetics and toxicity parameters of compounds is of paramount importance in reducing the late-stage attrition in drug development process. Top five compounds showed moderate solubility in water at 25 °C. High oral bioavailability can be expected from the top hits due to their high human intestinal absorption (⪆0.90%) [33]. All the hits also showed higher values for VDss than its upper limit (0.45) indicating these compounds may readily distribute into the body tissues rather than plasma. This broad tissue distribution is beneficial for antivirals and antibiotics [28]. Drugs that are either non-substrates or inhibitors of CYP450 enzymes can lead to accumulation of other drugs which pass via the same metabolic route. Therefore, understanding the drug metabolism is helpful for prevention of adverse drug-drug interactions which cause drug toxicity in treatments [34]. All the top hits used in this study were predicted to act as substrates of CYP3A4 but inhibitors for CYP1A2, CYP2C9, and CYP2C19. However, out of several CYP450 isoforms (CYP1A2/CYP2C9/CYP2C19/CYP3A4/CYP2D6) participating in detoxification of xenobiotics, CYP2D6 and CYP3A4 have the potential to metabolize 50% and 25% of known drugs, respectively [35]. The in silico results revealed that the selected top hits can undergo proper metabolism through the action of CYP3A4 enzyme. All the top hits also exhibited moderate total clearance (⪅ 0.51) allowing compounds to retain in the body for more extended periods at higher concentrations. Toxicity results of koenine and girinimbine revealed no harmful effect with regard to hepatotoxicity, carcinogenicity, mutagenicity, and cytotoxicity. Nevertheless, koenigicine, mukonicine, and o-methylmurrayamine A were shown to have carcinogenic and mutagenic potencies which can be modified through structure derivatization [36]. Overall, the top hits showed promising pharmacokinetic predictions along with less toxicity properties.

Previous research has reported on the pharmacological activities of the top hit phytochemicals studied here. Koenigicine was recognized to have antibacterial, antihyperglycemic, and antipancreatic lipase activities [37], whereas mukonicine was discovered to have antioxidant activity [10]. Girinimbine was also demonstrated to have antifungal, anti-inflammatory, anticancer, and antioxidant properties in several in vitro experiments [10]. In vitro antibacterial activities of koenine have been investigated [40]. Several in vitro tests have been conducted to study the pharmacological effects of o-methylmurrayamine A, including neuroprotective, antioxidant, and anti-inflammatory activities [10].

To date, no experimental results have been reported for the antiviral activity of the top hit compounds derived from Murraya koenigii. Against HCV and HIV, all of the top hits exhibited greater in silico antiviral inhibition percentages than 3WL. As a result, these hits are likely to be antiviral against SC2-Mpro. In support of this, recent in vitro studies demonstrated that some of clinically approved HIV and HCV drugs have potency to inhibit SC2-Mpro [38, 39].

Moreover, individual top hit phytochemicals and the top scoring pharmacophore model after alignment were inspected separately to identify the pharmacophoric features. In general, these pharmacophoric features are chemical features that are essential to interact with a specific biological target structure either to trigger or block its biological response. The examined top scoring pharmacophore model had shared structural features of 3 aromatic rings, 4 hydrophobic areas, 1 hydrogen bond donor, 1 hydrogen bond acceptor, and 1 positive ionizable region. Furthermore, this top-scoring pharmacophore model can be applied in ligand-based pharmacophore modeling in order to discover potent inhibitors against SC2-Mpro through independent database screening.

Conclusion

The docking findings of AutoDock 4.2 indicated five phytochemicals (koenigicine, mukonicine, o-methylmurrayamine A, koenine, and girinimbine) with greater binding affinities and lower inhibition constants than the reference inhibitor (3WL). They also demonstrate a distinct binding style and stronger interactions with the critical amino acid residues in the active region of the Mpro protein. Moreover, all top five phytochemicals showed strong binding to SC-2 Mpro of Alpha, Beta, Gamma, and Omicron variants. The top five phytochemicals are within the permissible range of physicochemical parameters and show satisfactory pharmacokinetic attributes. Each top hit molecule and best pharmacophore model has a different set of pharmacophoric features: aromatic ring, hydrophobic area, hydrogen bond donor/acceptor, and positively ionizable region. The MD simulation trajectories of RMSD, RMSF, and Rg of top five complexes exhibited high stability except koenigicine. Therefore, koenine, o-methylmurrayamine A, mukonicine, and girinimbine from Murraya koenigii, like most bioactive natural compounds, offer potential hits that may be further structurally modified and evaluated in vitro and in vivo for the discovery of novel SC-2 Mpro inhibitors.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 50 KB)

Supplementary file2 (XLSX 30 KB)

Supplementary file3 (DOCX 13 KB)