Computational prediction of nimbanal as potential antagonist of respiratory syndrome coronavirus.

The high pathogenic nature of the Middle East Respiratory coronavirus (MER) and the associated high fatality rate demands an urgent attention from researchers. Because there is currently no approved drug for the management of the disease, research efforts have been intensified towards the discovery of a potent drug for the treatment of the disease. Papain Like protease (PLpro) is one of the key proteins involved in the viral replication. We therefore docked forty-six compounds already characterized from Azadirachta indica, Xylopia aethipica and Allium cepa against MERS-CoV-PLpro. The molecular docking analysis was performed with AutoDock 1.5.6 and compounds which exhibit more negative free energy of binding, and low inhibition constant (Ki) with the protein (MERS-CoV-PLpro) were considered potent. The physicochemical and pharmacokinetic properties of the compounds were predicted using the Swissadme web server. Twenty-two of the compounds showed inhibition potential similar to dexamethasone and remdesvir, which had binding affinity of -6.8 and -6.3 kcal/mol respectively. The binding affinity of the compounds ranged between -3.4 kcal/mol and -7.7 kcal/mol whereas; hydroxychloroquine had a binding affinity of -4.5 kcal/mol. Among all the compounds, nimbanal and verbenone showed drug likeliness, they did not violate the Lipinski rule neither were they inhibitors of drug-metabolizing enzymes. Both nimbanal and verbenone were further post-scored with MM/GBSA and the binding free energy of nimbanal (-25.51 kcal/mol) was comparable to that of dexamethasone (-25.46 kcal/mol). The RMSD, RMSF, torsional angle, and other analysis following simulation further substantiate the efficacy of nimbanal as an effective drug candidate. In conclusion, our study showed that nimbanal is a more promising therapeutic agent and could be a lead for the discovery of a new drug that may be useful in the management of severe respiratory coronavirus syndrome.

Introduction

The Middle East respiratory syndrome coronavirus (MERS-CoV) was discovered to be highly pathogenic with potentials to infect human [ [12,45]] after it was detected in a Saudi Arabia man in 2012 [8,45]. Due to international travels of infected people, MERS-CoV has spread worldwide with 2502 laboratory-confirmed infection cases reported between September 2012 and December 2019 as well as 858 associated deaths [30]. In December 2019 a novel coronavirus-severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) formerly known as the 2019 novel coronavirus (2019-nCoV) and similar to SARS-CoV was identified in Wuhan, China [19,40]. Although SARS-CoV, MERS-CoV, and SARS-CoV-2 were identified to be highly pathogenic in the human population, there is presently no effective remedy against the virus [30].

The high case-fatality rate (CFR) of MERS-CoV infection far exceeds that of all other known human coronaviruses, including the human severe acute respiratory syndrome coronavirus (SARS-CoV) [21], therefore making it of great concern.

MERS-CoV similar to SARS-CoV is a single-strand positive-sense RNA virus whose large polyprotein is processed by two proteases, i.e. a 3-C-like protease (3CLpro) and a papain-like protease (PLpro) [21]. The single MERS-CoV papain-like protease [18,37] enzyme is part of a large non-structural protein 3 (nsp3) that contain four other domain; a ubiquitin like fold (UB1), an ADP-ribose-1d-phosphatase (ADRP) domain, a SARS-unique domain (SUD), and a transmembrane (TM) domain. PLpro enzyme function in the cleavage of the first three positions of its polyprotein, while 3CLpro cleaves the remaining 11 locations, releasing 16 nonstructural proteins (nsp) in both MERS-CoV and SARS-CoV [21]. Aside from the cleavage function, the MERS-PLpro also exhibits deubiquitination and de-ISGylation [25,43]. This deubiquitination and de-ISGylation function makes PLpro an interferon antagonist. It blocks the interferon regulatory factor 3 (IRF3) pathway [25,44], and evades the host cells innate immune response [10,25,43].

The use of plants for the management of various kinds of diseases constitutes great interest globally [39]. Many drug-like molecules present in plants are still of value in providing support to human health worldwide [32]. The current attention given to plant derived anti-microbial agents is linked to their safety and elongated history of practice [35]. Azadirachta indica, Xylopia aethipica and Allium cepa have been reported with various biological functions such as antiviral, antibacterial, antimalarial, anticancer and anti-inflammatory [3,4,9,11,20,36,42]. Therefore, they could serve as a repository of important biological compounds with inhibitory potentials against respiratory syndrome coronavirus. In this study, molecular modelling techniques that involve the combination of molecular docking, molecular dynamics simulations, MM/GBSA computation, and pharmacokinetics studies were employed for the characterization of target compounds already identified in these plants. These methods put together are robust and will provide a lead for the identification of promising drug candidates.

Methods

Protein preparation

The three-dimensional structure of MERS-CoV Papain like protease (PLpro) (PDB: 4PT5), was retrieved from the Protein Data Bank (PDB) (https://www.rcsb.org/) and prepared for molecular docking using BIOVIA Discovery Studio 2019. Before docking, the protein was prepared by removing interacting ligands, and water molecules downloaded together with it followed by saving the clean protein in the PDB format.

Ligand preparation

We obtained compounds already identified in Allium cepa, Azadirachta indica or Xylopia aethiopica through literature search ( Table 4 ) and extracted the SDS format of these compounds from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). We took the extracted SDS into the cactus online smiles translator (https://cactus.nci.nih.gov/translate/) for 3D PDB ligand download.

Table 4

Classification of compounds according to the plant source and their structure.

Compound source	Compound name	Chemical ID	Compound structure
A.indica	Meliacinin	15885442	Image 2
	Nimbanal	14194023	Image 3
	Nimbionol	189704	Image 4
	Nimbionone	189706	Image 5
	Nimbolide	100017	Image 6
	Nimocinol	178770	Image 7
	Quercetin	5280343	Image 8
	Regorafenib	11167602	Image 9
	Salannol acetate	14194026	Image 10
	3-deacetylsalanin	14458886	Image 11
	Azadirachtin	5281303	Image 12
	Azadironic acid	15885443	Image 13
A. cepa
	Cycloallin	12305353	Image 14
	Gamma-s-propyl-cysteine	13598411	Image 15
	Isohamnetin 3,4-diglucoside	5901757	Image 16
	Isorhamnetin 4-glucoside	44259381	Image 17
	Isorhamnetin	5281654	Image 18
	Isovallinin	12127	Image 19
	Kaempferol-3-O-rutinoside	5318767	Image 20
	Luteolin	5280445	Image 21
	Apigenin	5280443	Image 22
	Methiin	9578071	Image 23
	Quercetin 3,4-diglucoside	5320835	Image 24
	Quercetin 3,7,4-triglucoside	44259184	Image 25
	Quercetin 3-glucoside	5280804	Image 26
	Quercetin 7,4-diglucoside	11968881	Image 27
	Quercetin	5280343	Image 28
	Rutin	5280805	Image 29
	s-propylcysteine	125198	Image 30
	Zwiebelane	29322215	Image 31
	Apigenin	5280443	Image 32
X. aethiopica
	Bornyl acetate	6448	Image 33
	Carvone	7439	Image 34
	Citral	638011	Image 35
	Citronellol	8842	Image 36
	Copaene	12303902	Image 37
	Cryptone	92780	Image 38
	Cubebene	91747196	Image 39
	Cuminal	326	Image 40
	Methy chavicol	8815	Image 41
	Myrtenal	61130	Image 42
	Terpiene-4-ol	11230	Image 43
	Thymol	6989	Image 44
	Trans carveol	94221	Image 45
	Verbenone	29025	Image 46
	Alpha terpineol	17100	Image 47

Amino acid sequence alignment and analysis

We retrieved the FASTA sequence of seven viral proteins from the NCBI portal on June 30, 2020 and took the retrieved sequence for multiple sequence alignments. The list of the protein accession numbers is in Table 1 . BLASTP was used for the protein sequence alignments and the Phylogenetic tree was constructed using the Neighbor-Joining Method in Mega-X.

Table 1

Accession and Description of seven MERS PLpro.

Accession	Description	Links
Image 1 4PT5_A	Chain A, Papain-like Protease [Human betacoronavirus 2c EMC/2012]	Related Information
Image 1 4RNA_A	Chain A, papain-like protease [Human betacoronavirus 2c EMC/2012]	Related InformationStructure-3D structure displaysIdentical Proteins-Identical proteins to 4PT5_A
Image 1 5V6A_A	Chain A, MERS-CoV PLpro [Human betacoronavirus 2c EMC/2012]	Related InformationStructure-3D structure displaysIdentical Proteins-Identical proteins to 4REZ_A
Image 1 5V69_A	Chain A, MERS-CoV PLpro [Human betacoronavirus 2c EMC/2012]	Related InformationStructure-3D structure displaysIdentical Proteins-Identical proteins to 4REZ_A
Image 1 4RF1_A	Chain A, Orf1ab Protein [Human betacoronavirus 2c Jordan-N3/2012]	Related InformationStructure-3D structure displaysIdentical Proteins-Identical proteins to 4REZ_A
Image 1 4RF0_A	Chain A, Orf1ab Protein [Human betacoronavirus 2c Jordan-N3/2012]	Related InformationStructure-3D structure displaysIdentical Proteins-Identical proteins to 4REZ_A
Image 1 4REZ_A	Chain A, Orf1ab Protein [Human betacoronavirus 2c Jordan-N3/2012]	Related InformationStructure-3D structure displaysIdentical Proteins-Identical proteins to 4REZ_A

Molecular docking

The selected compounds were docked against MERS-CoV-PLpro (PDB ID: 4PT5) to examine the molecular interactions existing between the docked complex. Blind docking and docking calculation were achieved using AutoDock Tool 1.5.6 and AutoDock Vina respectively [38]. We achieved protein clean-up and visualization of molecular interaction between receptor and ligand with BIOVIA Discovery studio 2019. Formation of stable complex between the protein and ligand indicates the high potency of an inhibitor [17,29]. Compounds, which exhibit more negative free energy of binding and low Ki, calculated using equation (1) were more potent. The hydrogen bonds and hydrophobic interactions between the ligands and the protein were also studied using BIOVIA Discovery Studio 2019. Identity of compounds obtained from each plant species is revealed in Table 4 . During docking, we added polar-H-atoms to the protein followed by Gasteiger charges calculation. The protein file was saved as pdbqt file and the grid dimensions were set as 84, 66, 100 in x, y and z directions while the Centre grid box was set as −6.584, −5.089, and 3.027 for x, y and z respectively.where ΔG is the binding energy in kcal/mol, R is the universal gas constant (8.314 J/mol/K) and T is the temperature (298 K).

ADME profiling

Drug likeliness and Pharmacokinetics prediction of the compounds was achieved computationally using the Swissadme web server ( www.swissadme.ch/ ). The SMILES structure of the screened compounds obtained from the PubChem database was used for the ADME prediction.

Binding free energy calculation/thermodynamics calculation

We applied the molecular mechanics/Poisson-Boltzmann surface area (MM/GBSA) module integrated with the prime of the Schrodinger suite for the calculation of binding free energy of docked complex output and computed the relative free energy of the ligand complex using the OPLS3 force field, VSGB solvent and the rotamer search algorithm [28].

The binding free energy was extrapolated from the equation below Eqn (2);

Molecular dynamics simulation studies

Desmond module of Schrodinger software was employed for the molecular dynamics simulation [31]. We selected nimbanal and dexamethasone because of their good binding affinity and pharmacokinetic properties for the MD simulation with OPLS 2005 force field parameters. The docked complexes were centred on the orthorhombic box of the predefined TIP3P water system. The box's volume was minimized, and the net charge of the system was neutralized by incorporating 0.15 M NaCl into each system to model physiological state [5]. The temperature and pressure were kept constant at 300 K and 1.01325 bar using Nose–Hoover thermostat and Martyna–Tobias–Klein barostat methods. Simulation analysis was performed through the NPT ensembles by considering heavy atoms, time intervals, and pressure [24,46]. 20 ns simulation of the relaxed system were carried out with NPT ensembles, and the long–range electrostatic interactions were computed using the Particle–Mesh–Ewald algorithm. The trajectories were recorded at 4.8 ps intervals, and the protein-ligand interaction, stability, and behaviour was performed using the Desmond simulation interaction diagram in maestro [1,31].

Results and discussion

Sequence alignments

Of the seven proteins, three were of Jordan while the other four were of England origin (Fig. 1 ). The sequence alignments of 4PT5 template protein with the six other MERS-CoV protein only showed 100% identity with the PDB structure 4RNA_A, whereas absolute identity were observed among the protein structures with the accession numbers 4REZ_A, 4RF0_A, 4RF1_A, 5V69_A and 5V6A_A. The alignment analysis showed that all the proteins have the same number of amino acid residue and that a mutation might have occurred at the position 100 and 113 in 4PT5 and 4RNA to warrant variation from 4REZ_A, 4RF0_A, 4RF1_A, 5V69_A and 5V6A_A (Fig. 2 ).

Fig. 1

Phylogenic tree of seven Middle East Respiratory Syndrome Coronavvirus Papain Like Protease.

Fig. 2

Multiple Sequence Alignment of six MERS-CoV with the template 4PT5.

Molecular docking study

For the validation of the docking protocol, the root mean square deviation obtained for the ligands before and after docking was close to 1, therefore the docking protocol is correct.

The Middle East Respiratory Syndrome coronavirus Papain like protease consist of 324 amino acid sequence with molecular weight of 36.10 kDa. Previous study stated that the protein is made of N-terminal ubiquitin-like (UBI) domain and the catalytic domain (Lei et al., 2014). The catalytic domain comprises of three distinct subdomains i.e. the thumb, fingers and palm. The substrate-binding site lies between the palm and thumb domains whereas the catalytic traid consisting of Cys 111, His 278 and Asp 239 is at the centre (Lei et al., 2014).

Residues Leu106–Tyr112, Gly161–Arg168 in the thumb subdomain, Phe269–Tyr279, Pro250, and Thr308 in the palm subdomain (Fig. 4 ) of the MERS-CoV PLpro, line the substrate-binding site. Replacement of Tyr 279, corresponding to Tyr274 in SARS-CoV PL pro enzyme by Ala leads to a loss of protease activity (Barretto et al., 2005; Lei et al., 2014). The X-ray diffraction structure of MERS-CoV-PLpro, PDB:4PT5 has a resolution of 2.59 Å, Ramachandran outlier of 1.3% and Rfree value of 0.298.

Fig. 4

(A) Cartoon view of the MERS-CoV PLpro enzyme's overall structure. α-Helices (cyan) and β-strands (purple) are numbered, polypeptide segments devoid of repetitive secondary structure, including loops and turns, are brown. The ubiquitin-like (Ubl) domain is encircled by a red dashed line. The catalytic domain consists of the thumb, fingers, and palm subdomains. A gray sphere indicates the structural zinc ion in the fingers domain. The Cα atoms of the catalytic-site cysteine (111), histidine (278), and aspartate (293) residues are also shown (yellow, blue, and red sphere, respectively). The red arrow indicates the substrate-binding region and points to the catalytic site. (B) The four-cysteine ligands (Cys191, Cys194, C226 and C228) and the structural zinc ion (gray sphere) in the zinc ribbon of the fingers domain. Sulfur atoms are shown in yellow, oxygen in red, nitrogen in blue, and carbon in light blue. (C) The catalytic triad: Cys111, His278, and Asp293. (Adapted from Ref. [22]. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

The selected compounds from Allium cepa, Azadirachta indica or Xylopia aethiopica were docked against the Middle East Respiratory Syndrome coronavirus (MERS-CoV-PLpro (PDB:4PT5) (Fig. 3 ) in order to examine the non-bonding interaction present in the studied complex. The potency of the natural compounds was benchmarked on the binding affinity of the standard inhibitors; hydroxychloroquine, dexamethasone and remdesvir. The binding affinity of the compounds ranged from −3.4 to −7.7 kcal/mol. Only 22 of the compounds were shortlisted when their binding affinities and inhibition constant were compared with the standards. The shortlisted compounds either compare favourably with remdesvir (−6.3 kcal/mol) or dexamethasone (−6.8 kcal/mol), i.e. compounds with binding affinity ≤ −6.0 kcal/mol, were excluded. Remdesivir is a broad-spectrum antiviral drug currently under clinical trial. It has demonstrated activity against RNA viruses in several families, including Coronaviridae (such as SARSCoV, MERS-CoV), and strains of bat coronaviruses capable of infecting human respiratory epithelial cells [7,33,34]. The drug was reported to have some inhibitory effects against coronavirus replication [41].

Fig. 3

3D-View of Middle East Respiratory Syndrome coronavirus.

Quercetin 3-glucoside was the best among the compounds requiring the lowest energy to interact with PLpro (inhibition constant (Ki) of 2.7 μM). The order of antiviral potency of the compounds includes; Quercetin 3-glucoside > Quercetin 7,4-diglucoside > Azadironic acid = Kaempferol-3-O-rutinside = Quercetin 3,7,4-triglucoside > Quercetin = Quercetin 3,4-diglucoside > Isorhamnetin 4-glucoside = Isohamnetin 3,4-diglucoside > Regorafenib = Quercetin-4-glucoside > Rutin = Nimbolide = Luteolin = Nimocinol > Isorhamnetin > Meliacinin = Apigenin > Nimbionone > Nimbanal > Cubebene > Azadirachtin. Quercetin 3-glucoside, Quercetin 7,4-diglucoside, Kaempferol-3-O-rutinside, Quercetin 3,7,4-triglucoside, Quercetin, Quercetin 3,4-diglucoside, Isorhamnetin 4-glucoside, Isohamnetin 3,4-diglucoside, Quercetin-4-glucoside. Rutin, Luteolin, Isorhamnetin, Apigenin were identified in Allium cepa while Azadironic acid, Meliacinin, Nimbionone, Nimbanal, Azadirachtin, Regorafenib have been identified in Azadirachta indica and Cubebene in Xylopia aethiopica (Table 4).

Quercetin-3-glucoside showed stronger affinity for PLpro than others, having binding affinity of −7.7 kcal/mol and forming hydrogen bond interaction with ASP 65, ASN 16, VAL 10 and ASP 13 and hydrophobic interaction with ALA 69 and LYS 68. Nimbolide with −6.9 kcal/mol as binding affinity and 10.4 μm as inhibition constant formed hydrogen bond interaction with THR 67, LEU 62, PRO 79, ALA 59 and LEU 82 but hydrophobic contact with LEU 70.

Nimbanal only formed hydrogen bond interaction with one residue, LYS 143 of MERS-CoV-PLpro, but formed hydrophobic constant with LYS 287, LEU 124, ASP 123, LYS 126, HIS 142, AGR 285, LYS 102, VAL 99, VAL 103 and ARG 104. Each of the compounds were either involved in one or more of hydrophobic and hydrogen bond interaction with the amino acid residue of the protein (Table 2). The docking pose of some of the top ranked docked compounds is displayed in Fig. 5 A to N. The strong negative calculated binding affinity of 22 of the natural compounds indicates they have promising inhibitory effect and could serve as starting point in the development of effective drugs targeting pulmonary respiratory viral diseases.

Table 2

Docking score and the interactions of the natural compounds with PLpro.

S/N	Ligands	Binding Affinity ΔG (Kcal/mol)	Inhibition Constant Ki (μM) 10−6	Interacting Amino acids	Bond Type
	Hydroxychloroquine	−4.5	562.5	VAL 217, GLU 218, ALA 222, TYR 319, ARG 223, SER 321, ASP 322, ASN 324, GLN 237, ARG 221	Van der waals, Conventional Hydrogen Bond, Pi-Pi T-shaped, Pi-Alkyl
	Remdesvir	−6.3	28.2	ASP 13, LEU 73, VAL 15, ALA 69. LYS 68, GLU 72, ASP 75, PHE 17, VAL 10, ARG 18, THR 11, THR 63, GLU 66, ASN 16	Van der waals, Attractive Charge, Conventional Hydrogen bond, Carbon hydrogen Bond, Pi-Cation, Pi-Anion, Pi-Sigma, Alkyl, Pi-Alkyl
	Dexamethaxone	−6.8	12.3	ASN 159, VAL 77, TYR 157, ILE 132, LYS 178, PRO 76, HIS 173, VAL 212	Conventional Hydrogen bond, unfavorable bump, Alkyl, Van der waals,
	3-deacetylsalanin	−6.2	33.3	VAL 212, TYR 211, LYS 207, ALA 177, LEU 176, LEU 206, HIS 173, CYS 210, MET 187, LEU 203, ARG 236,	Van der waals, conventional hydrogen bond, unfavorable donor-donor, alkyl, Pi-alkyl
	Alpha terpineol	−4.4	664.3	GLU 72, PRO 133, VAL 15, LEU 73, ASN 16, ALA 69, ASP 13	Van der waals, Alkyl
	Apigenin	−6.4	23.9	LYS 71, LEU 70, THR 67, LEU 62, ASN 61, ASP 60, ALA 59, PRO 79, LEU 82, ASP 78, VAL 77	Van der waals, conventional hydrogen bond, Carbon hydrogen bond, Pi-Sigma, Pi-Alkyl
	Azadirachtin	−6.0	46.4	VAL 99, MET 97, GLY 145, ARG 104, TRP 95, ASP 146, THR 148, ILE 151, SER 147, LYS 96, HIS 93, GLY 94	Van der waals, conventional hydrogen bond, Pi-alkyl
	Azadironic acid	−7.3	5.35	PRO 79, PRO 76, LYS 71, THR 67, PHE 81, LEU 70, ASP, 78, VAL 77, LEU 82, ALA 59, LEU 62, ASP 60, ASN 61	Van der waals, conventional hydrogen bond, Alkyl
	Bornyl acetate	−5.0	244.9	HIS 142, LYS 143, ARG 104, ASP 123, ARG 285, MET 120, VAL 103, LYS 102	Van der waals, conventional hydrogen bond, Alkyl, Pi-Alkyl
	Buoebenone	−5.7	76.5	GLU 218, THR 261, THR 258, VAL 257, PRO 263, LYS 255, PHE 265, PHE 314, THR 216, VAL 217, ASP 264	Van der waals, Alkyl, Pi-Alkyl
	Carvone	−5.2	175.6	LYS 126, PRO 244, PRO 315, TRP 245, PHE 314, LEU 313, VAL 266, LEU 124	Van der waals, Alkyl
	Citral	−5.4	125.9	THR 216, PHE 268, GLN 215, VAL 309, VAL 312, PHE 265, VAL 257, THR 258, THR 261, PRO 263, LEU 256, LYS 255	Van der waals, conventional hydrogen bond, Alkyl, Pi-Alkyl
	Citronellol	−4.7	403.4	THR 67, ALA 59, LEU 62, ASN 61, ASP 60, PRO 79, ASP 78, PHE 81, LEU 82, VAL 77, LEU 70	Van der waals, conventional hydrogen bond, Alkyl, Pi-Alkyl
	Copaene	−5.6	90.3	PHE 81, ASP 78, VAL 77, LEU 70, LEU 82, THR 67, PRO 79, LEU 62, ASP 60, ALA 59	Van der waals, Alkyl
	Cryptone	−4.5	562.5	LEU 70, ASP 78, PHE 81, LEU 82, ALA 59, PRO 79, LEU 62, VAL 77	Van der waals, conventional hydrogen bond, Alkyl
	Cubebene	−6.0	46.4	THR 261, ASP 264, GLU 218, PHE 265, VAL 257, THR 258, LYS 255, PRO 263, THR 216, LEU 256	Van der waals, lkyl, Pi-Alkyl
	Cuminal	−5.3	148.7	ASN 61, LEU 62, PRO 79, VAL 77, LEU 70, ASP, 78, PHE 81, THR 80, LEU 82, ALA 59	Van der waals, conventional hydrogen bond, Alkyl, Pi-Alkyl
	Cycloallin	−4.8	341.6	PHE 85, LYS 89, ASP 149, LYS 141, ALA 138, HIS 137, LEU 73, ALA 134, ASP 13, TYR 74, TYR 57, TYR 86,	Van der waals, conventional hydrogen bond, Alkyl, Unfavorable donor-donor
	Gamma-s-propyl-cysteine	−4.0	1291.9	LEU286, VAL 266, GLY 289, THR 259, GLY 288, LEU 290, ALA 262, PRO 263, ASP 264, PRO 315, PHE 265	Van der waals, conventional hydrogen bond, Carbon hydrogen bond, Alkyl
	Isohamnetin 3,4-diglucoside	−7.1	7.5	LYS 141, ASP 149, ASP 146, THR 148, LYS 89, TYR 57, ASP 13	Conventional hydrogen bond, carbon hydrogen bond, Pi-cation, Pi-Pi Stacked, Pi-Pi T-shaped
	Isorhamnetin 4-glucoside	−7.1	7.5	ILE 128, PHE 139, MET 140, LYS 129, CYS 182, ASP 127, LEU 125, LYS 126, LYS 143, HIS 142, ASP 123, ARG 104, MET 97, ARG 285, VAL 103, LYS 103	Van der waals, conventional hydrogen bond, carbon hydrogen bond, unfavorable donor-donor, Pi-Cation, Pi-Ananion, Pi-Pi Stacked, Pi-Alkyl
	Isorhamnetin	−6.5	20.2	ALA 59, ASP 60, LEU 82, THR 67, PRO 76, LYS 71, LEU 62, LEU 70, VAL 77, PHE 81, ASP 78, PRO 79	Van der waals, conventional hydrogen bond, unfavorable acceptor-acceptor, Pi-Sigma, Pi-Alkyl
	Isovallinin	−4.7	403.4	LEU 70, LEU 62, PRO 79, LEU 82, THR 80, PHE 81, ASP 78, VAL 77, THR 67	Van der waals, conventional hydrogen bond, Pi-Sigma, Pi-Alkyl
	Kaempferol-3-O-rutinside	−7.3	5.3	ALA 134, LEU 73, TYR 74, VAL 12, ASP 13, GLY 14, TYR 57, HIS 93, THR 148, ASP 146, LYS 141, LYS 89, ASP 149, ALA 90, GLY 39, PHE 36, PHE 37, PHE 85, TYR 86, ASN 38	Van der waals, conventional hydrogen bond, Pi-Cation, Pi-Pi T-shaped, Pi-Alkyl
	Luteolin	−6.9	10.4	ASN 61, THR 67, VAL 77, PHE 81, ASP 78, HIS 83, LEU 82, THR 80, ALA 59, PRO 79, LEU 70, LEU 62	Van der waals, conventional hydrogen bond, Unfavorable donor-donor, Pi-Alkyl
	Meliacinin	−6.4	23.9	PHE 265, LEU 286, LYS 287, PHE 314, LYS 126, TRP 245, PRO 244, LEU 124, LEU 313, PRO 315, VAL 266	Van der waals, conventional hydrogen bond, Alkyl, Pi-Alkyl
	Methiin	−3.7	2127.5	PHE 81, ASP 78, THR 80, LEU 82, HIS 83, PRO 79, ALA 59, LEU 62, THR 67, VAL 77, LEU 70	Van der waals, conventional hydrogen bonds
	Methy chavicol	−5.1	207.4	ALA 59, ASP 60, LEU 62, ASN 61, THR 67, LEU 82, PRO 79, LEU 70, VAL 77, ASP 78, PHE 81	Van der waals, conventional hydrogen bond, Alkyl, Pi-Sigma
	Myrtenal	−5.2	175.6	TYR 74, ASP 149, PHE 85, TYR 86, LYS 89, LEU 153, LEU 73, ALA 138, ALA 134, HIS 137, LYS 141	Van der waals, conventional hydrogen bond, Alkyl, Pi-Alkyl
	Nimbanal	−6.0	46.4	VAL 99, ARG 104, LYS 143, HIS 142, LYS 126, ASP 123, LEU 124, LYS 287, ARG 285, VAL 103, LYS 102	Van der waals, Carbon hydrogen bond, Pi-Cation, Alkyl, Pi-Alkyl
	Nimbionol	−5.7	76.5	LYS 102, VAL 99, ARG 285, ASP 123, ARG 104	conventional hydrogen bond, Carbon hydrogen bond, Alkyl, Pi-Alkyl
	Nimbionone	−6.2	33.4	THR 67, ASN 61, LEU 62, ALA 59, PRO 79, LEU 82, ASP 78, PHE 81, VAL 77, LEU 70	Van der waals, Carbon hydrogen bond, Unfavorable Acceptor-Acceptor, Alkyl, Pi-Alkyl
	Nimbolide	−6.9	10.4	ASP 78, LEU 82, ALA 59, PRO 79, LEU 70, LEU 62, THR 67	conventional hydrogen bond, Carbon hydrogen bond, Alkyl, Pi-Alkyl
	Nimocinol	−6.8	12.3	LEU 176, LEU 206, CYS 210, HIS 173, LYS 207, VAL 212, TYR 211, MET 187, LEU 203, TRP 189, ARG 236	Van der waals, Pi-Sigma, Alkyl, Pi-Alkyl
	Quercetin 3,4-diglucoside	−7.2	6.3	TYR 57, LYS 141, ASP 149, LYS 89, ASN 38, HIS 93, SER 55	Conventional hydrogen bond, Pi-Pi T-shaped
	Quercetin 3,7,4-triglucoside	−7.3	5.3	LYS 129, CYS 182, ILE 128, ASP 127, LYS 143, LEU 125, LEU 122, LYS 126, PHE 139, ARG 139, VAL 103, LYS 102, ASP 123, ARG 104, LEU 124, HIS 142, TRP 245, VAL 266, LEU 313, LEU 246, PHE314, PRO 315, GLY 316	Van der waals, conventional hydrogen bond, Carbon hydrogen bond, Unfavorable donor-donor,Pi-Alkyl
	Quercetin 3-glucoside	−7.7	2.7	LYS 68, ALA 69, ASP 65, ARG 18, VAL 18, ASN 16, ASP 13	conventional hydrogen bond, Carbon hydrogen bond, Pi-Sigma, Amide-Pi Stacked Pi-Alkyl
	Quercetin 7,4-diglucoside	−7.6	3.2	HIS 173, LYS 207, GLU 180, SER 184, CYS 183, ALA 185, GLY 205, GLN 204, LEU 181, LEU 206, ALA 177, LEU 176, CYS 210	Van der waals, conventional hydrogen bond, Carbon hydrogen bond, Pi-Cation, Pi-Sigma, Amide-Pi Stacked, Pi-Alkyl
	Quercetin	−7.2	6.3	PHE 81, VAL 77, ASP 78, LEU 70, PRO 79, LEU 82, LEU 62, THR 67, ASN 61, ALA 59, ASP 60	Van der waals, conventional hydrogen bond, Unfavorable donor-donor, Pi-Alkyl
	Quercetin-4-glucoside	−7.0	8.8	SER 184, ARG 186, ALA 185, CYS 183, GLN 204, GLY 205, LEU 206, LYS 207, LEU 176, ALA 177, ALA 179, GLU 180, LEU 181	Van der waals, conventional hydrogen bond, Carbon hydrogen bond, Pi-Cation, Pi-Alkyl
	Regorafenib	−7.0	8.8	GLU 180, LYS 207, GLY 205, LEU 181, ALA 179, LEU 206, LEU 176, TYR 211, VAL 212, HIS 173, CYC 210, ALA 177	Van der waals, conventional hydrogen bond, Halogen (Fluorine), Pi-Anion, Pi-Sulfur, Pi-Alkyl
	Rutin	−6.9	10.4	ASN 61, ASP 60, ALA 59, PRO 79, THR 80, PRO 79, LEU 82, ASP 78, PHE 81, LEU 70, VAL 77, GLY 75, PRO 76, LYS 71, LEU 62, THR 67	Vam der waals, conventional hydrogen bond, Unfavorable acceptor-acceptor, Pi-Anion
	Salannol acetate	−5.5	106.6	ALA 177, HIS 173, VAL 212, ARG 236, LYS 207, CYS 210	Conventional hydrogen bond, Alkyl, Pi-Alkyl
	s-propylcysteine	−3.4	3503.8	MET 120, ASP 123, ARG 104, HIS 142, LYS 102, VAL 103, ARG 285	Van der waals, conventional hydrogen bond, Alkyl
	Terpiene-4-ol	−4.4	664.3	LEU 73, PRO 133, VAL 15, HIS 137, ASN 16, ALA 69, ASP 13, GLU 72	Van der waals, conventional hydrogen bond, Alkyl, Pi-Alkyl
	Thymol	−5.4	125.9	ALA 59, VAL 77, LEU 70, LEU 82, PHE 81, ASP 78, PRO 79, LEU 62, THR 67, ASP 60	Van der waals, conventional hydrogen bond, Alkyl, Pi-Alkyl
	Trans carveol	−5.1	207.4	ASP 13, PHE 85, ALA 134, HIS 137, LYS 141	conventional hydrogen bond, Alkyl, Pi-Alkyl
	Verbenone	−5.3	148.7	LYS 141, ALA 138, HIS 137, ALA 134, LEU 73, TYR 74	conventional hydrogen bond, Alkyl, Pi-Alkyl
	Zwiebelane	−4.3	784.4	TYR 74, ALA 134, LEU 73, ASP 13, HIS 137, LYS 141, ALA 138, ASP 149, PHE 85, LEU 153	Van der waals, Alkyl, Pi-Alkyl

Fig. 5

A–N: 3D and 2D complex structure of binding between ligand and protein.

ADME prediction study

The drug-attrition rate in preclinical and clinical trials is quite alarming due to poor pharmacokinetic studies [15]. Consequently, initial screening of these drug-like molecules can increase the chances of passing through the clinics [16]. In this regard, the docked compounds were subjected to in silico ADME screening. The drug-likeness test was based on the Lipinski's Rule of Five [23]. The criteria for the Lipinski's violation was premise on the distribution of the compound, molecular weights (MW), calculated lipophilicity (log P), number of hydrogen bond acceptors (HBA) and number of hydrogen bond donors (HBD); (MW < 500; log P < 5; HBD≤5; HBA≤10) ([23,27]. Apigenin, Bornyl acetate, Isorhamnetin, Luteolin, Nimbionone, Nimbionol, Nimbolide, Nimocinol, Quercetin, Regorafenib Verbenone and Hydroxychloroquine have zero Linpinski violation (Table 3b). In addition, these compounds together with 3-deacetylsalanin and Azadironic acid are highly absorbed by the intestine whereas only Verbenone, Nimbionone, Nimbionol, Nimbolide, Cubebene and Buoebenone, Bornyl acetate, Copaene and Hydroxychloroquine have blood-brain barrier permeation. The cytochrome P450 enzymes (CYPs) constitute a superfamily of proteins important in the metabolism and detoxification of xenobiotics [6]. Inhibition of any of the drug-metabolizing CYPs will elevate the concentration of the corresponding drug substrate and bring about drug overdose [26]. Some of the compounds with probable drug likeliness including hydroxychloroquine and remdesivir were inhibitors to one or more of the metabolizing enzymes Table 3a . Interestingly, Nimbanal, Nimbolide, Nimocinol and Verbenone were non-inhibitor to any of the metabolizing enzymes, therefore these compounds are potential drug with attractive pharmacokinetic profiles. P-glycoprotein (Pgp) is a critical determinant of the pharmacokinetic properties of drugs as it functions to extracts foreign substances from the cell [2]. Nimbanal and Verbenone were predicted to be non P-glycoprotein (Pgp) substrate and therefore appeared to be the only among all to have drug-like potentials.

Table 3b

Lipinski violation on molecular weights (MW), calculated lipophilicity (log P), number of hydrogen bond acceptors (HBA) and number of hydrogen bond donors (HBD).

Compounds	MW	HBA	HBD	Molar refractivity	Mlog P	Lipinski Violations
Hydroxychloroquine	335.87	3	2	98.57	2.35	0
Remdesvir	602.58	12	4	150.43	0.18	2
3-deacetylsalanin	554.67	8	1	147.07	2.54	1
Apigenin	270.24	5	3	73.99	0.52	0
Azadironic acid	454.6	5	1	129.05	4.21	1
Azadirachtin	720.71	16	3	165.92	−0.47	2
Cubebene	204.35	0	0	67.14	5.65	1
Buoebenone	204.35	0	0	67.14	5.65	1
Bornyl acetate	196.29	2	0	56.33	2.76	0
Copaene	204.35	0	0	67.14	5.65	1
Isohamnetin 3,4-diglucoside	640.54	17	10	146.75	−4.42	3
Isorhamnetin 4-glucoside	478.4	12	7	114.63	−0.24	2
Isorhamnetin	316.26	7	4	82.5	−0.31	0
Kaempferol-3-O-rutinside	594.52	15	9	139.36	−3.43	3
Luteolin	286.24	6	4	76.01	−0.03	0
Meliacinin	512.68	6	0	144.11	4.38	2
Nimbanal	510.58	8	0	132.91	1.81	1
Nimbionone	302.26	4	1	84.17	1.89	0
Nimbionol	304.38	4	2	85.13	1.98	0
Nimbolide	466.52	7	0	120.99	2.28	0
Quercetin 3,4-diglucoside	626.52	17	11	142.28	−2.7	3
Quercetin 3,7,4-triglucoside	788.66	22	14	174.4	−6.64	3
Nimocinol	452.58	5	1	126.44	3.36	0
Quercetin 3-glucoside	464.38	12	8	110.16	−2.59	2
Quercetin	302.24	7	5	78.03	−0.56	0
Quercetin 7,4-diglucoside	626.52	17	11	142.28	−4.62	3
Regorafenib	482.82	8	3	112.44	3.28	0
Rutin	610.52	16	10	141.38	−3.89	3
Salannol acetate	598.72	9	0	157.28	2.94	1
Verbenone	150.22	1	0	45.42	2.2	0

Table 3a

In silico evaluation of the ADME profile of the selected compounds.

Compounds	GI Absorption	BBB Permeant	P-gp Substrate	CYPIA2 Inhibitor	CY2C19 Inhibitor	CYP2C9 Inhibitor	CYP2D6 Inhibitor	CYP3A4 Inhibitor
Hydroxychloroquine	High	Yes	No	Yes	No	No	Yes	No
Remdesvir	Low	No	Yes	No	No	No	No	Yes
3-deacetylsalanin	High	No	Yes	No	No	No	No	No
Apigenin	High	No	No	Yes	No	No	Yes	Yes
Azadironic acid	High	No	Yes	No	No	Yes	No	Yes
Azadirachtin	Low	No	Yes	No	No	No	No	No
Cubebene	Low	Yes	No	Yes	Yes	Yes	No	No
Buoebenone	Low	Yes	No	Yes	Yes	Yes	No	No
Bornyl acetate	High	Yes	No	No	No	Yes	No	No
Copaene	Low	Yes	No	Yes	Yes	Yes	No	No
Isohamnetin 3,4-diglucoside	Low	No	Yes	No	No	No	No	No
Isorhamnetin 4-glucoside	Low	No	Yes	No	N0	No	No	No
Isorhamnetin	High	No	No	Yes	No	No	Yes	Yes
Kaempferol-3-O-rutinside	Low	No	Yes	No	No	No	No	No
Luteolin	High	No	No	Yes	No	No	Yes	Yes
Meliacinin	High	No	Yes	No	No	Yes	No	No
Nimbanal	High	No	No	No	No	No	No	No
Nimbionone	High	Yes	Yes	No	No	No	No	Yes
Nimbionol	High	Yes	Yes	No	No	No	Yes	No
Nimbolide	High	No	Yes	No	No	No	No	No
Quercetin 3,4-diglucoside	Low	No	Yes	No	No	No	No	No
Quercetin 3,7,4-triglucoside	Low	No	Yes	No	No	No	No	No
Nimocinol	High	No	Yes	No	No	No	No	No
Quercetin 3-glucoside	Low	No	No	No	No	No	No	No
Quercetin	High	No	No	Yes	No	No	Yes	Yes
Quercetin 7,4-diglucoside	Low	No	Yes	No	No	No	No	No
Regorafenib	Low	No	No	Yes	Yes	Yes	Yes	Yes
Rutin	Low	No	Yes	No	No	No	no	No
Salannol acetate	Low	No	Yes	No	No	No	No	Yes
Verbenone	High	Yes	No	No	No	No	No	No

MM/GBSA binding free energy

Measurement of binding strength and affinity of a ligand when it occupies the active region of the protein using thermodynamics calculation will help to understand the nature of existing interaction [13]. Both verbenone and nimbanal showed interesting pharmacokinetic properties and together with dexamethasone, were post-scored with MM/GBSA. This method of scoring is well established to show a more reliable statistical relationship to experimental binding affinity [14]. From the MM/GBSA computation, the binding free energy for dexamethasone, nimbanal and verbenone was −25.46, −25.51 and −9.14 kcal/mol respectively (Fig. 10). The result further justified that the inhibitory potential of nimbanal was comparable to that of dexamethasone and so, only the two compounds were subjected to molecular dynamics simulation.

Fig. 10

Docking score (Binding Affinity) and the MM/GBSA Binding energy of dexamethasone, nimbanal and verbenone.

Molecular dynamics simulations of dexamethasone and nimbanal with MERS-CoV-PLpro

Nimbanal, among the natural compounds, exhibited drug like characteristics and strong binding affinity, and was selected with dexamethasone for MD simulation. The dynamic simulation was performed for 20 ns, and the dynamics stability of the complexes was analysed using parameters like the L-RMSF, protein-ligand RMSD, protein-ligand contacts, ligand torsion profile, P-RMSF and ligand properties.

Root mean square deviations of protein (RMSD-P), ligand root mean square fluctuation (L-RMSF) and protein root mean square Fluctuation (P-RMSF)

To determine the conformational stability of the protein backbone and the protein-ligand complex, the RMSD, which measures the distance between the protein backbones of the superimposed protein, was monitored. With the P-RMSD, we evaluated when simulation has equilibrated by studying the movements of different atoms in the protein when in contact with the ligand at the active site. From Fig. 6 a, the RMSD-P for dexamethasone-protein complex ranged between 1.5 and 4.0 Å while the value ranged from 1.0 to 4.0 Å when in complex with nimbanal. (Fig. 6b). This observation showed that nimbanal better maintained the stability of the protein throughout the simulation period. To assess how stable the ligand was to the protein-binding pocket, the Lig fit Prot was observed and it shows the RMSD of the ligand when the protein-ligand was first aligned on the protein backbone. When the L-RMSD values are larger than the RMSD of the protein, it is an indication that the ligand has diffused away from its initial binding site. The L-RMSD value of dexamethasone was ranged between 0.5 and 3.0 Å (Figs. 6a) and 1.8–14 Å for nimbanal (Fig. 6b). The Lig fit on Prot (RMSD-P and RMSD-L) showed that the MERS-CoV-PLpro ligand complex is stable over the molecular dynamics simulation time.

Fig. 6

Plot of RMSD values of Dexamethasone and Nimbanal (a, b respectively) with Mpro RMSD value as a function of time. The protein RMSD is at the Left Y-axis and the ligand RMSD is at the right Y-axis coordinates. L-RMSF for dexamethasone and nimbanal (c, d respectively) and P-RMSF of dexamethasone and nimbanal (e, f respectively).

The RMSF plot further indicates the areas of the protein that fluctuate the most during the MD simulation. Customarily, the tails (N- and C- terminal) fluctuate the most in the protein structure. Conversely, the Secondary Structure Elements (SSE) like the alpha helices and beta strands are usually more rigid, therefore, fluctuate less than the loop regions. The MERS-CoV-PLpro was made of 19.41% alpha helix, 32.41% beta strand to make 51.82% SSE (Fig. 9). The RMSF of the MERS-CoV-PLpro ligand complex (L-RMSF) and the MERS-CoV-PLpro (P-RMSF) are shown in Fig. 6c–f. Dexamethasone atoms from 1 to 3 highly fluctuated (RMSF >0.75 Å, Fig. 6c) whereas, nimbanal atoms ranging from 26 to 33 (RMSF>6 Å, Fig. 6d) highly fluctuated. In Fig. 6e and f, each peak indicates the protein area that fluctuates the most during the course of MD simulation. For all the amino acids in MERS-CoV-PLpro when in complex with dexamethasone and nimbanal, the P-RMSF were below 4.8 Å Fig. 6e and f.

Fig. 9

Protein secondary structure elements (SSE) (red-alpha helices, blue-beta strands) of Mpro with the reports of SSE distribution by residue and the summary of the SSE composition of each trajectory frame over the course of simulation. In the bottom each residue and its SSE assignment over time is shown. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Protein-ligand interactions (or 'contacts') are categorized into four types: Hydrogen Bonds, Hydrophobic, Ionic, and Water Bridges. In Fig. 7 a–e, interaction that existed between the atoms of ligand and protein residues more than 30.0% of the simulation time in the selected trajectory (0.00 through 20.02 nsec) are shown. A Hydrophobic bond was observed between TYR 211, VAL 212 and hydroxyl (OH) on C 28 of dexamethasone. OH, group of C 27 also formed a polar bond with HIS 173 and O group of C 2 formed charged (positive) interaction with ARG 236. These interactions were observed to be stable over 98%, 69%, and 31%, respectively, of the simulation time (Fig. 7a). Nimbanal on the other hand formed water bridges with the protein. In addition, each interaction type contains specific subtypes that can be explored through the ‘Simulation Interactions Diagram’ panel (Fig. 7e and f). Normalization of the stacked bar charts was observed over the course of the trajectory: For example, a value of 0.7 suggests that in 70% of the simulation time, a specific interaction is maintained. Values over 1.0 are possible as some protein residue may make multiple contacts of the same subtype with the ligand. In the bar diagram, the fraction of interactions with each amino acid residue over the course of the simulation run is specified (Fig. 7c and d).

Fig. 7

A schematic of detailed dexamethasone and nimbanal (a, and b respectively) atoms interactions with the amino acid residues'. Bar charts of protein interaction with dexamethasone and nimbanal (c and d respectively) as monitored throughout the simulation (green-H-bonding; gray-hydrophobic; blue-water bridges; pink-ionic interactions). Plot of the contacts and interactions between protein and ligand: dexamethasone and nimbanal (e and f respectively) over the course of trajectory. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

The ligand torsions plot summarizes the conformational evolution of every rotatable bond (RB) in the ligand throughout the simulation trajectory. The top panel shows the 2D schematic of a ligand with colour-coded rotatable bonds. Each rotatable bond torsion is accompanied by a dial plot and barplots of the same colour. Dial (or radial) plots describe the conformation of the torsion throughout the course of the simulation. The beginning of the simulation is in the centre of the radial plot and the time evolution is plotted radially outwards. The bar plots summarize the data on the dial plots, by showing the probability density of the torsion. If torsional potential information is available, the plot also shows the potential of the rotatable bond (by summing the potential of the related torsions). The values of the potential are on the left Y-axis of the chart and are expressed in kcal/mol. Finally, looking at the histogram and torsion potential relationships may give insights into the conformational strain the ligand undergoes to maintain a protein-bound conformation (Fig. 8 a–d ).

Fig. 8

2D schematic plot of dexamethasone and nimbanal with colour-coded rotatable bond (a, and b respectively). The dial plots and bar plots of rotatable bond torsions for dexamethasone and nimbanal (c and d respectively). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Conclusion

From this study, nimbanal and verbenone appear to be the only compounds with attractive physicochemical and pharmacokinetic properties to qualify as drug-like molecule. In addition, nimbanal appeared to be a good inhibitor for MERS-CoV-PLpro with high binding affinity and low inhibition constant. From our observation, the high binding affinity of nimbanal correlates with it stable interaction with the protein following simulation. However, non nimbanal-protein interaction exhibited might suggest further in-vivo or in-vitro analysis to establish the pharmacological and biological properties of nimbanal. Overall, the simulation result showed that nimbanal is a more promising therapeutic agent in the treatment of severe respiratory syndrome compared to dexamethasone. The favourable interaction of dexamethasone might contribute to its importance in the pool of library of trial drugs. Furthermore, the RMSD, RMSF, torsional angle, and other analysis following simulation demonstrated that nimbanal could be an effective drug candidate. Therefore, it provides a potential lead for the treatment of severe respiratory syndrome.

Funding

(The authors received no funding from any source).

Availability of data and material

(Not applicable).

Code availability

(Not applicable).

Authors' contributions

(API, SB and AAE conceptualize the study, API, FOS, AA, AZA provided the resources, API and AAE performed the docking study, AAE did the ADMET study, ATS and API performed the BLAST analysis, API, AA, FOS, AZA and BTA did the molecular dynamics study, API and BTA did the MM/GBSA study, API wrote the manuscript, all authors approved the final draft).

Ethics approval

(Not applicable).

Consent to participate

(Not applicable).

Consent for publication

(All authors approved the manuscript for publication).

Declaration of competing interest

The authors categorically declare that there exist no financial or any conflict of interest whatsoever among them.