Molecular docking, simulation and MM-PBSA studies of nigella sativa compounds: a computational quest to identify potential natural antiviral for COVID-19 treatment.

Nigella sativa or black seed is used as a medicinal plant around the globe. Oil and seeds have a long tradition of folklore use in various medicinal and food systems. The conventional therapeutic use of Nigella sativa, in different ways, has been reported in several studies to treat different diseases including influenza, headache, hypertension, diabetes, inflammation, eczema, fever, cough, asthma, bronchitis, and fever. Based on previously reported potential therapeutic uses of N. sativa compounds, and keeping in mind the dire need of time for the development of potent antiviral, a combined docking, ADMET properties calculation, molecular dynamics, and MM-PBSA approaches were applied in the current study to check the therapeutic potentials of N. sativa chief constituents against COVID-19. Among the studied compounds, we found that dithymoquinone (DTQ), with binding affinity of -8.6 kcal/mol compared to a positive control (chloroquine, -7.2 kcal/mol) , has the high potential of binding at SARS-CoV-2:ACE2 interface and thus could be predicted as a plausible inhibitor to disrupt viral-host interactions. Molecular dynamics simulation of 100 ns well complemented binding affinity of the compound and revealed strong stability of DTQ at the docked site. Additionally, MM-PBSA also affirms the docking results. Compound DTQ of the present study, if validated in wet lab experiments, could be used to treat COVID-19 and could serve as a lead in the future for development of more effective natural antivirals against COVID-19. Communicated by Ramaswamy H. Sarma.

Introduction

Nigella sativa, a herbal medicine belonging to the Ranunculaceae family, has been explored exclusively and acquired worldwide acknowledgment as a proven historical and religious remedy for a wide range of health problems (Goreja, 2003). N. sativa seeds have been consumed to treat various diseases and infirmities (Begum & Mannan, 2020; Majeed et al., 2020; Yimer et al., 2019). In Islamic as well as Christian literature it is considered as one of the best curing medicine. In “Tibb-e-Nabwi” the regular use of black seed is recommended. Unlike antibacterial drugs, the availability of antiviral medicines is limited with narrow effectiveness (D’Alessandro et al., 2020). Oil extracted from black seed was used effectively in the “Murine cytomegalovirus (MCMV) viral inhibition where reduced titer in the liver and spleen of infected individual was reported (Salem & Hossain, 2000). Black seed has been used in turkey poultry to significantly suppress the pathogenicity of the influenza virus (H9N2) and also resulted in improved immune response (Umar et al., 2016). N. sativa seeds showed strong antiviral properties when used against “Laryngotrachietis Virus (ILTV)” at a concentration of 35 μM in Chicken Embryo Rough Cells (CER) (Zaher et al., 2008). In a clinical trial conducted on Egyptian patients who were infected with the Hepatitis C virus, some of them were diabetic and HCV-positive and were also not recommended for IFN/ribavirin therapy, N. sativa was given to them on regular basis for nearly three months (450 mg thrice a day). On completion of treatment, an improvement in the clinical condition of patients with decreased viral road, less oxidative stress and improved glycemic control was observed (Barakat et al., 2013). Additionally, another surprising study was reported in 2013, in which an adult who was HIV-positive was treated with N. sativa at a dose of 10 ml (twice a day) for the 6 months. After 6 months, when the “EIA Western Bolt” test was repeated, the outcome was seronegative, the test was repeated multiple times, and results were obtained all the time, affirming the efficacy and potency of N. sativa seeds as potent anti-HIV (Onifade et al., 2013). The medical importance of N. stativa has led to increased research interest to test its therapeutic potential against the COVID-19.

The rapid outbreak of COVID-19 in Wuhan china followed by its spread worldwide with nearly 1,016,395 infected people and thousands of deaths, poses a serious public health issue (Boopathi et al., 2020; Elflein, 2020). COVID-19 is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) with the mild symptoms such as dry cough and rise in body temperature to severe symptoms like pneumonia or acute respiratory failure, particularly in adults (males) with comorbidity (Chen et al., 2020; Zhou et al., 2020). Coronavirus most likely originated from animals probably bats, and transmission occurs to humans that progressed in the form of infected droplets (Chen et al., 2020; Wu et al., 2020; Zhou et al., 2020). COVID-19 has been declared as a global health emergency by the World Health Organization (WHO), and it is believed that this pandemic like SARS-CoV-1, would result in the loss of many lives with tremendous social effects and huge economic damages (Keogh-Brown & Smith, 2008). No specific and potent therapeutics are reported against this pandemic yet (Aanouz et al., 2020; Chen et al., 2020; Khan et al., 2020). Thus, the designing and development of effective antivirals along with vaccines against SARS-CoV-2 is the dire need of the world today (Bouchentouf & Missoum, 2020; Emeka et al., 2014; Enayatkhani et al., 2020; Khan et al., 2020; 2020; Muralidharan et al., 2020; Sarma et al., 2020). Keeping the current situation in mind and diverse therapeutically use of N. sativa, the present investigation aimed to identify the active constituents of N. sativa as probable drug leads against COVID-19.

Material and methods

Ligands preparation

The 2D structures of N. sativa chief constituents were generated using a freely available ChemDraw software of Chem Office 2004 (Li et al., 2004). The structures were then converted into .mol2 format with the help of Chem3D Ultra, another package of Chem Office. For structural geometry optimization of selected inhibitor freely available UCSF Chimera 1.14 was employed (Pettersen et al., 2004).

Target selection and preparation

Recently reported three dimensional (3D) structure of the COVID-19 chimeric receptor-binding domain with its human receptor ACE2 (PDB ID: 6VW1) was retrieved from RCSB PDB (https://www.rcsb.org/structure/6vw1). The water molecules as well as co-crystallized ligands were deleted from the PDB file.

Docking protocol

To determine the binding mode and interaction of the selected compounds and target, docking studies were performed using Autodock/vina (Trott & Olson, 2010). The pdbqt files of the receptor protein, and N. sativa compounds along with the grid box setting at the active site of the receptor for compounds binding was done through Auto Dock GUI program (Morris et al., 2009). The grid size boundaries along X, Y, and Z axes was scaled at 40 Å with a grid spacing of 0.375 Å to allow proper binding flexibility at the docked site. The output pdbqt files were written into a configuration (conf) file. The receptor was treated as rigid entity whereas ligands were kept flexible to attain the best fitting conformation with respect to the receptor complex. The generated solutions of docking were clustered and those with root mean square deviation (RMSD) value <1.0 Å were considered only. The binding conformation of ligands with the lowest binding affinity was characterized as the most stable conformation of the ligands with respect to the receptor. The docking protocol was validated by docking previously crystalized co-ligands with the S-protein. The visualization tools including Visual Molecular Dynamics (VMD) (Humphrey et al., 1996), and Chimera were employed to analyze intermolecular chemical interactions such as both hydrogen bonding, ionic interactions, and hydrophobic contacts.

ADME properties analysis

The ADME properties of the selected compounds were predicted by freely available online SwissADME software (Daina et al., 2017). Toxicity analysis was performed through online ProTox-II software (Drwal et al., 2014). Toxicity parameters, such as hepatotoxicity, carcinogenicity, immunogenicity, mutagenicity, and cytotoxicity were evaluated. Toxicity is important to evaluate in drug designing as it helps in determining the toxic dose in animal model studies and also lessens the number of animal model studies (Sripriya et al., 2019).

Molecular dynamic (MD) simulation

The selected complex then underwent 100 ns MD simulation using AMBER 18 software (Case et al., 2018). The Leap of Amber19 tools was used to create receptor molecule topology files (Case et al., 2005). The simulation system was neutralized by adding Na + ions. The neutralized system was then submerged in TIP3PBOX, a water molecules box (Jorgensen et al., 1983). The solvated system was minimized first for 1500 steepest descent method iterations, followed by 1000 steps of conjugate gradient (Andleeb et al., 2016). System heating was achieved for 100 ps that involve gradual heating from 0 K to 300 K at pressure of 1 atm (Abbasi et al., 2016). After, the system was equilibrated for 100 ps (Lavenda, 2016) that allowed the exchange between kinetic and potential energies. Simulation production run was accomplished for 100 ns to evaluate dynamics of the complex and check stability of ligand docked conformation. To apply constraint on covalently bound hydrogen atoms, SHAKE algorithm (Ryckaert et al., 1977) was employed. The periodic boundary conditions in the simulation box were ensured with the use of canonical ensemble (Ahmad et al., 2019). Temperature of the was kept constant using Berendsen algorithm (Lemak & Balabaev, 1994). MD simulations were completed using the Ewald summation method (Darden et al., 1993) and the coordinate files were saved for every 0.5 ps. The trajectories were analyzed using the CPPTRAJ module implemented in Amber18.

Binding free energies of the complex

The free energy of binding for the docked complex was estimated using two methods: Molecular Mechanics-Generalized Born Surface Area (MMGBSA) and Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) embedded in MMPBSA.py module of AMBER18 (Miller et al., 2012). In total, 100 frames were processed from the trajectories and the net energy of the system was estimated through the following equation,

Each of the above terms involves the calculation of several energy components: van der Waals energy, electrostatic energy, and internal energy summed from molecular mechanics and polar contribution towards solvation energy. The analysis also takes into account the contribution from non-polar term towards solvation energy and inhibitor entropy.

Results

Nigella sativa and derived compounds

Active ingredients of N. sativa reported in several studies are dithymoquinone (DTQ), thymoquinone (TQ), thymol, thymohydroquinone (THQ), p-cymene, 4-terpineol, and t-anethole. N. sativa seeds also contain other ingredients such as vitamins, mineral elements, carbohydrates, fats, proteins, and essential amino acids (Ahmad et al., 2013; Al-Jassir, 1992; Bhatia & Bajaj, 1972; Chun et al., 2002; Correa et al., 1986; Islam et al., 2019; Majeed et al., 2020; Omar et al., 1999). Bioactive compounds used in the current study are shown in Figure 1.

Figure 1.

Two-dimensional structures of compounds studied: (a) 4-TERPINEOL (b) Dithymoquinone (DTQ) (c) Carvacrol (d) Carvone (e) p-cymene (f) t-anethole (g) Thymohydroquinone (THQ) (h) Thymol (THY) and (i) Thymoquinone (TQ),.

Molecular docking

Docking studies were conducted triplicate to have a clear depiction of the mode of action of N. sativa chief constituents against S protein: ACE2 receptor interface. For each compound, ten interactions were generated and the one with best binding affinity was selected. Docking results revealed that among the studied compounds, DTQ exhibited a high affinity of −8.6 kcal/mol for the selected target (Table 1).

Table 1.

Binding affinities of docked compounds using Autodock/vina.

Compounds	Binding Affinities (kcal/mol)	Distance from Best Mode
		RMSD Lower Bound	RMSD Upper Bound
A-Terpineol	−5.8	0.000	0.000
Dithymoquinone	−8.6	0.000	0.000
Carvacrol	−7.0	0.000	0.000
Carone	−6.5	0.000	0.000
P-Cymene	−5.8	0.000	0.000
T-Anethole	−6.2	0.000	0.000
Thymoquinone	−6.7	0.000	0.000
Thyhydromoquinone	−6.1	0.000	0.000
Thymol	−6.0	0.000	0.000
Chloroquine (positive control)	−7.2	0.000	0.000

Interacting residues of the interface of S protein RBD (PHE 28, LYS31, GLU35, GLN76, LEU79) and ACE2 receptor complex (LYS452, LEU455, PHE456, GLU484, GLY485, CYC488, TYR489, PHE490, LEU492, GLN 493and SER494), involved in ligand binding, are shown in Figure 2. Detailed 3 D inhibitor- interface interactions are depicted in Figure 3, which may be useful for future identification of novel inhibitors and for designing of analogs of the studied compound. Exploration of the docked complex using VMD revealed that DTQ formed three hydrogen bonds and one ionic interaction with a key interacting residue of interface LYS31 (LYS31:HZ3—UNK:O 1.80 Å, LYS31:HZ2—UNK:O 2.23 Å, LYS31:HZ1—UNK:O 2.68 Å, LYS31:NZ—UNK:O 1.79 Å). Ionic as well as hydrogen bond interactions were also observed with GLN493 (GLN493:NE2—UNK: O 3.201 Å, GLN493:2HE2—UNK:O 3.53 Å, GLN493:1HE2—UNK:O 3.16 Å), a residue important for binding with human ACE2. SER494 was also involved in electrostatics and H-bond interactions (SER494: HN—UNK:O 2.78 Å, SER494:N—UNK:O 3.60 Å). Hydrogen bond interactions play a vital role in maintaining the stability of the docked complex (Panigrahi, 2008). Additionally, a network of hydrophobic contacts between carbonyl groups of DTQ and interacting residues of target: LYS31:CD—UNK:C 1.18 Å, LYS31:CD—UNK:C 1.81 Å, LYS31:CD—UNK:C 2.41 Å, LYS31:CD—UNK:C 3.37 Å, LYS31:CD—UNK:C 3.66 Å, LYS31:CD—UNK:C 3.10 Å, LYS31:CG—UNK:C 1.62 Å, LYS31:CG—UNK:C 2.82 Å, LYS31:CG—UNK:C 3.84 Å, LYS31:CB—UNK:C 3.09 Å, LYS31:CB—UNK:C 3.49 Å, GLN76:CD—UNK:C 3.66 Å, TYR489:CE1—UNK:C 3.91 Å, TYR489:CD1—UNK:C 3.36 Å, TYR489:CB—UNK:C 3.94 Å, TYR489:CA—UNK:C 3.85 Å, PHE490:CD2—UNK:C 3.59 Å, PHE490:CD2—UNK:C 3.68 Å, PHE490:CB—UNK:C 3.67 Å, PHE490:CB—UNK:C 3.97 Å, PHE490:CG—UNK:C 3.92 Å, LYS452:CD—UNK:C 3.99 Å (Figure. 3). Aforesaid analysis of docked complex revealed that interactions noted between the DTQ and target complex have more preferences towards the ACE2 receptor, might in-turn limit the binding of S-protein with the human host, ACE2 receptor, and thus may restrict the viral interaction within the host cell. In addition, to further confirm docking results, we have superimposed DTQ-target docked complex with chloroquine-S-protein:ACE2 docked complex. Chloroquine (CQ) was recently tested in vitro and was found to be active against COVID-19. It was observed that more or less CQ and DTQ occupied the same binding pocket (Figure 4).

Figure 2.

Highlighting the key interacting residues of DTQ and S-protein: ACE2 active site generated using Chimera.

Figure 3.

Three-dimensional representation of important interactions using VMD.

Figure 4.

Superimposed structures of CQ and DTQ docked to S-protein:ACE2 interface.

To calculate the drug-likeness, the ADME properties of DTQ were calculated (Table 2). Based on the calculated LogP value, the compound proved to be lipophilic with a consensuses LogP value of 2.92. ESOL value (−3.20) indicates that the compound is moderately water soluble. DQT estimated to have decent absorption in the gastrointestinal tract making it a good drug candidate for oral route administration. No violation of Lipinski’s rule of five was noted (Lipinski, 2004). The feasibility of the suggested compound to serve as a lead for the future design of potent antivirals to treat COVID-19 was explored using medicinal chemistry filters. Toxicity analysis indicated DTQ as inactive for all the used toxicity prediction parameters and thus predicted as non-toxic.

Table 2.

ADME properties calculated for dithymoquinone using the SwissADME server.

ADME Properties	Values
Molecular weight	326.39 g/mol
Lipophilicity (WLogP)	3.09
Consensus Log Po/w	2.92
Water Solubility (ESOL)	−3.20
GI absorption	High
BBB permeant	Yes
Log Kp (skin permeation)	−6.41 cm/s
Lipinski rule	Yes, 0 violation
Leadlikeness	Yes
Synthetic accessibility	3.69

MD Simulation

The binding affinity and dynamics of top-ranked DTQ conformation obtained as a result of docking were further evaluated using MD simulation. RMSD was determined to check how the binding of DTQ to the active site of protein affects its ability to reach an equilibrium state. As shown in Figure 5, the system equilibrated after 0.5 ns, however, fluctuations were observed throughout the production run with an average of 2.58 Å. It was noticed that during the 100 ns simulation run, the ligand re-equilibrated nine times. The docked complex was snapshot at different intervals from the production trajectories to assess the behavior and stability of ligand, shown in Figure 6. The studied system experienced varying degrees of fluctuations at different intervals, i.e. around 12, 20, 29 (Figure 6a), 50, 52, 59 (Figure 6b), 76, 78, and 92 ns (Figure 6c). The superimposition of protein structures, extracted from the trajectories, revealed that at aforementioned intervals ligand moved a bit away from the binding site and after some time (about 0.5 ns) returned to its original position. Further, the stability of the docked complex was assessed by computing the mobility of residues. Analysis of the root mean square fluctuation (RMSF) graph revealed no major fluctuations in the active site region. The average estimated value was 1.16 Å. Lower RMSF values for the interacting residues indicated the loss of flexibility because of ligand binding as shown in Figure 7 (Ahmad et al., 2017; Kuzmanic & Zagrovic, 2010). Results have shown that these residues might play a vital role in stabilizing the complex.

Figure 5.

RMSD profile of DTQ and S-protein: ACE2 complex.

Figure 6.

Superimposition of structures extracted from trajectories at, (a) 12, 20, 29 ns, (b) 50, 52, 59 ns, (c) 76, 78, 92 ns interval of 100 ns MD simulation. Movements are shown by kaki, cyan, and purple, respectively.

Figure 7.

RMSF profile of DTQ and S-protein: ACE2 complex. Interacting residues are labeled red.

MMGBSA and MMPBSA binding free energies

The binding free energy of the complex was carried out to revalidate the inhibitor affinity for the receptor protein complex predicted by the docking simulation studies. The binding free energies for the system in both MMGBSA and MMPBSA are tabulated in Table 3. As can be seen, the complex interactions are mainly dominated and favored by van der Waals (−32.8290 kcal/mol) in both methods in contrast to less contribution from electrostatic (−4.7312 kcal/mol). The polar solvation energy was revealed more unfavorable in receptor-inhibitor formation (14.3225 kcal/mol in MMGBSA and 17.0286 kcal/mol in MMPBSA). The net binding energy for the complex predicted by MMGBSA is −26.7955 kcal/mol that revealed stronger system stability compared to the one concluded by MMPBSA (−22.9620 kcal/mol). The difference in the net energy of both techniques is large because of the polar solvation energy which seems to play less role in receptor-inhibitor interaction. The entropy impact over the inhibitor has not been considered due to higher computational cost.

Table 3.

Binding free energies for the DTQ and S-protein: ACE2 complex.

GENERALIZED BORN
Energy Component	Average	Std. Dev.	Std. Err. of Mean
VDWAALS	−32.8290	2.0424	0.2042
EEL	−4.7312	2.0910	0.2091
EGB	14.3225	1.7571	0.1757
ESURF	−3.5578	0.1562	0.0156
DELTA G gas	−37.5602	2.7240	0.2724
DELTA G solv	10.7647	1.7745	0.1774
DELTA TOTAL	−26.7955	1.8539	0.1854
POISSON BOLTZMANN
VDWAALS	−32.8290	2.0424	0.2042
EEL	−4.7312	2.0910	0.2091
EPB	17.0286	1.9171	0.1917
ENPOLAR	−2.4304	0.1289	0.0129
EDISPER	0.0000	0.0000	0.0000
DELTA G gas	−37.5602	2.7240	0.2724
DELTA G solv	14.5982	1.8847	0.1885
DELTA TOTAL	−22.9620	2.0078	0.2008

Discussion

An outbreak of COVID-19 in Wuhan in the Hubei province of China opens up new horizons for the exploration of synthetic as well as natural small compounds to identify a potential lead that could be used to design plausible antiviral drugs. Recently reported crystallographic structure of SARS-CoV-2 S-protein: human ACE2 receptor, indicates that disrupting the viral S-protein-ACE2 interface may be an attractive target for structure-based drug discovery (Elfiky, 2020a; 2020b; Enmozhi et al., 2020; Joshi et al., 2020). However, the discovery and development of new potent drugs is a time taking and costly process. It almost takes 10–15 years, including trials and regulatory approval, in the United States to report a novel drug (Van Norman, 2016). Keeping this and the reported therapeutic potentials of N. sativa main constituents in mind, the current in silico study is designed. Among the studies compounds, DTQ also known as nigellone, came out to be a more promising compound against the selected target. Nigellone, a carbonyl polymer of thymoquinone, retains its much of the pharmacologic properties (Chakravarty, 1993). Based on results of docking, dynamics, and net free energy value, we hypothesized that N. sativa compounds particularly DTQ (Eid et al., 2017; Shrivastava et al., 2011) is showing high affinity and stability at SARS-CoV-2:ACE2 interface which upon binding with key residues of the interface may disrupt host recognition and also cure the viral infection by disrupting the S-protein pathway. Beside SARS-CoV-2:ACE2 interface, several other proteins could also target for drug designing (Elfiky & Azzam, 2020; Pant et al., 2020). Among these, RNA-dependent RNA polymerase (RdRp) (catalyzes the replication of RNA from RNA template) (Elfiky, 2020b), ion channel E protein (Gupta et al., 2020), and the main protease (MPro) (Joshi et al., 2020) that is vital in viral replication and maturation can serve potential drug targets against SARS-CoV-2.

Conclusions

In conclusion, N. sativa chief constituents revealed to have potential as natural antiviral agents against novel coronavirus. Particularly, DTQ showed a high affinity towards SARS-CoV-2:ACE2 interface and interacts with several hotspot residues both through hydrophobic and hydrophilic bonding. Additionally, this compound has high solubility, gut absorption, and good profile of druglikness making it a potent lead for future structure optimization to design more potent derivatives. Dynamically in an aqueous environment, the compound predicted conformation of docking studies is highly stable at CoV-2:ACE2 interface as affirmed by the binding free energy assays by securing a very low energy level. Based on these findings, it is hypothesized that DTQ (PubChem SID: 163237941) may be subjected to in vitro and in vivo assays validation to disclose its true inhibitory potency.