In silico investigation of spice molecules as potent inhibitor of SARS-CoV-2.

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a novel infectious disease that is in rapid growth. Several trials are going on worldwide to find a solution for this pandemic. The viral replication can be blocked by inhibiting the receptor-binding domain (RBD) of SARS-CoV-2 spike protein (SARS-CoV-2 RBD Spro) and the SARS-CoV-2 main protease (SARS-CoV-2 Mpro). The binding of potential small molecules to these proteins can inhibit the replication and transcription of the virus. The spice molecules that are used in our food have antiviral, antifungal and antimicrobial properties. As spice molecules are consumed in the diet, hence its antiviral properties against SARS-CoV-2 will benefit in a significant manner. Therefore, in this work, the molecular docking of 30 selected spice molecules (screened through ADME property) was performed to identify the potential inhibitors for the RBD Spro and Mpro of SARS-CoV-2. We have found that though all the molecules bind actively with the SARS-CoV-2 RBD Spro and Mpro, but Piperine has the highest binding affinity among the 30 screened molecules. Besides, the comparative study between Piperine and currently used drugs show that Piperine is more effective. The interaction of Piperine with RBD Spro and Mpro is further validated by the molecular dynamics (MD) simulation studies. The free energy landscape and binding free energy results also, support for the stable complex formation of Piperine with RBD Spro and Mpro. We anticipate immediate wet-lab experiments and clinical trials in support of this computational study that might help to inhibit the SARS-CoV-2 virus. Communicated by Ramaswamy H. Sarma.

Introduction

The novel coronavirus disease 2019 (COVID-19) has become a major threat worldwide due to its fast-spreading nature. This disease is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The entry of this coronavirus to the host cell is mediated through the transmembrane spike glycoproteins (Hasan et al., 2020). This glycoprotein consists of two subunits and is reported to have a similar affinity to the human angiotensin-converting enzyme 2 (ACE2) as that of the severe acute respiratory syndrome coronavirus (SARS-CoV), which in turn results in efficient spreading of SARS-CoV-2 in humans (Walls et al., 2020). The spike glycol protein binds to its receptor human ACE2 by its receptor-binding domain (RBD) and is activated proteolytically by human protease (Shang et al., 2020). The interaction of the RBD of the spike glycoprotein to ACE2 is carried out by ARG403, TYR453, SER494, TYR495, PHE497, GLN498, THR500, ASN501, TYR505 residues of spike glycoprotein (Lan et al., 2020). The interaction of RBD Spro to ACE2 can be inhibited by the small molecules that interact with the above residues of RBD spike protein (Spro). On the other hand, the SARS-CoV-2 main protease (Mpro), also known as chymotrypsin-like protease, or 3-C-like protease (3CLpro), plays a vital role in processing the polyproteins through the translation of viral RNA. This protease is reported to have a minimum of 11 cleavage sites resulting in viral replication and toxicity (Zhang et al., 2020). The inhibition of these two viral targets can actively block the fusion and replication of SARS-CoV-2.

Currently, researchers are working globally for finding treatment for this disease in identifying a specific drug or vaccine that can inhibit viral replication at the earliest possible. The devastation of this disease is vividly seen from the data on the WHO website, which shows the infected patient number more than 245 lakhs and casualties more than 8 lakhs, worldwide from 216 countries, and still, it is continuing (https://www.who.int/emergencies/diseases/novel-coronavirus-2019, 29 August 2020). Currently, there are no approved drugs and vaccines for the treatment of COVID-19, but a few drugs such as remdesivir, hydroxychloroquine, etc. are under restricted use in case of emergency (Magagnoli et al., 2020). Meanwhile, the computational tools, molecular docking and molecular dynamics (MD) have gained attention as essential tools to investigate potential inhibitor molecules (Anurag et al., 2020; Gupta et al., 2020; Sourav et al., 2020). Choy et al. reported in-vitro studies showing remdesivir, lopinavir, emetine and homoharringtonine inhibits SARS-CoV-2 replication (Choy et al., 2020). Similarly, Wang et al. reported the inhibition property of remdesivir and chloroquine against novel coronavirus (Wang et al., 2020). In addition to different drug compounds, researchers also searched for natural molecules having antiviral activity. Natural constituents from foods, spices, herbs are also being found to have anti-infective properties. In this context, small active molecules present in natural products and their derivatives have gained tremendous attention as a source of therapeutic agents due to structural diversity for many years.

From 1940 to 2014, the US Food and Drug Administration (FDA) has approved about 49% of all small molecules that are natural products or derivates linked directly to those molecules (Newman & Cragg, 2016). The compound of essential garlic oil, a spice used in food, is reported as an inhibitor using the molecular docking method (Thuy et al., 2020). There are several recent studies on the inhibition of SARS-CoV-2 using many different natural and antiviral molecules (Al-Khafaji et al., 2020; Joshi et al., 2020; D. Kumar, Kumari, et al., 2020; Muralidharan et al., 2020). Recently, Das et al., using blind molecular docking, investigated for the potential inhibitors of SARS-CoV-2 Mpro (Sourav et al., 2020). Molecules studied by Das et al. are drug molecules, antivirals, antifungals, anti-nematodals and anti-protozoals in addition to natural compounds. Besides, natural molecules such as alkaloids and terpenoids from African medicinal plants were studied by Gyebi et al. for the inhibition property against SARS-CoV-2 Mpro (Gyebi et al., 2020). Recently, Umesh et al. screened compounds from Indian spices as potent inhibitors of SARS-CoV-2 Mpro (Umesh et al., 2020). Every spice has a particular aroma, colour and flavour due to the presence of specific molecules in them, and also, have antiviral properties (Aboubakr et al., 2016; Astani et al., 2010; Brochot et al., 2017; Chang et al., 2013; Choi, 2016; Mair et al., 2016; Zhang et al., 2014). These properties of the spice molecules compel us to conduct the present study, where we investigated the inhibition property of molecules present in various spices against the SARS-CoV-2 RBD Spro and SARS-Cov-2 Mpro using molecular docking and MD simulation studies. The compounds tested and their source of origin with PubChem ID are listed in .

Table 1.

Predicted data of docking score, solubility, pharmacokinetics, drug-likeness and medicinal chemistry of the screened spice molecules.

Sl. No.	Molecule	Binding energy (kcal/mol)		ESOL Log S	Ali Log S	Silicos-IT LogSw	GI absorption	BBB permeant	Pgp substrate	Log Kp (cm/s)	Bioavailability score	Synthetic accessibility
		Mpro	RBD Spro
1	2-Decenoic acid	–5.4	–4.6	–2.8	–4.23	–2.15	High	Yes	No	–4.68	0.56	2.44
2	α-Terpinyl acetate	–5.5	–5.0	–3.35	–4.21	–2.36	High	Yes	No	–4.69	0.55	3.13
3	Capsaicin	–6.4	–5.5	–3.53	–4.5	–4.87	High	Yes	No	–5.62	0.55	2.32
4	Carvone	–6.2	–5.2	–2.41	–2.72	–2.16	High	Yes	No	–5.29	0.55	3.33
5	Cinnamaldehyde	–5.7	–5.1	–2.17	–1.88	–2.4	High	Yes	No	–5.76	0.55	1.65
6	Cuminaldehyde	–5.9	–5.1	–2.52	–2.37	–3.15	High	Yes	No	–5.52	0.55	1.0
7	Dipropyl disulfide	–3.1	–3.0	–2.14	–3.42	–2.48	High	Yes	No	–5.3	0.55	2.79
8	Eucalyptol	–5.2	–4.9	–2.52	–2.59	–2.45	High	Yes	No	–5.3	0.55	3.65
9	Linalool	–5.5	–4.9	–2.4	–3.06	–1.84	High	Yes	No	–5.13	0.55	2.74
10	Vanillin	–5.7	–4.8	–1.82	–1.78	–1.88	High	Yes	No	–6.37	0.55	1.15
11	Thymol	–5.8	–5.3	–3.19	–3.4	–3.01	High	Yes	No	–4.87	0.55	1.0
12	Sabinene hydrate	–5.2	–4.7	–2.07	–2.18	–1.91	High	Yes	No	–5.74	0.55	2.82
13	Piperine	–7.3	–6.4	–3.74	–3.96	–3.0	High	Yes	No	–5.58	0.55	2.92
14	Menthol	–5.6	–5.2	–2.88	–3.5	–1.48	High	Yes	No	–4.84	0.55	2.63
15	Eugenol	–6.0	–5.0	–2.46	–2.53	–2.79	High	Yes	No	–5.69	0.55	1.58
16	Estragole	–5.7	–4.8	–3.09	–3.24	–3.35	High	Yes	No	–4.81	0.55	1.28
17	Gingerol	–6.1	–5.5	–2.96	–3.82	–4.58	High	Yes	No	–6.14	0.55	2.81
18	Shogaol	–5.8	–5.4	–3.7	–4.67	–4.8	High	Yes	No	–5.15	0.55	2.51
19	Paradol	–6.0	–4.6	–3.72	–4.79	–5.52	High	Yes	No	–5.08	0.55	2.28
20	Zingerone	–6.0	–5.1	–1.8	–1.68	–3.1	High	Yes	No	–6.7	0.55	1.52
21	Borneol	–5.7	–4.3	–2.51	–2.8	–1.91	High	Yes	No	–5.31	0.55	3.43
22	Bornyl acetate	–5.3	–4.8	–3.63	–4.57	–2.58	High	Yes	No	–4.44	0.55	3.64
23	Citral	–5.5	–4.7	–2.43	–3.05	–1.96	High	Yes	No	–5.08	0.55	2.49
24	Citronellal	–4.8	–4.8	–2.88	–3.88	–2.33	High	Yes	No	–4.52	0.55	2.57
25	2-Undecanone	–4.9	–4.3	–2.94	–4.15	–3.83	High	Yes	No	–4.43	0.55	1.72
26	Geranyl acetate	–5.4	–4.8	–3.21	–4.3	–2.52	High	Yes	No	–4.63	0.55	2.72
27	Nerolidol	–5.8	–5.0	–3.8	–4.99	–3.15	High	Yes	No	–4.23	0.55	3.53
28	Terpinen-4-ol	–5.2	–5.5	–2.78	–3.36	–1.91	High	Yes	No	–4.93	0.55	3.28
29	Terpineol	–5.7	–5.2	–2.87	–3.49	–1.69	High	Yes	No	–4.83	0.55	3.24
30	Decanal	–4.7	–3.9	–2.67	–3.85	–3.44	High	Yes	No	–4.56	0.55	1.62

Materials and methods

In general, molecular modelling is implemented as an essential tool for the prediction of drug–macromolecule interaction. This technique helps to enhance the success rate of an experiment and cuts down the experimental cost. Hence, the molecular docking study can help to analyse the possible binding pose of a small molecule on the active site of a macromolecule. Here we used molecular docking to screen some biologically active spice molecules with the SARS-CoV-2 RBD Spro and SARS-Cov-2 Mpro. The molecule with the highest binding affinity to RBD Spro and Mpro was subjected to MD simulation for further validation.

Drug-likeness properties of the small molecules

The property of the small molecules for drug-likeness was estimated using the Lipinski’s rule (Lipinski et al., 2001). This rule works on five parameters viz. no more than five hydrogen bond donors, no more than 10 hydrogen bond acceptors, molecular mass <500 Da, and the octanol-water partition coefficient, i.e. log P should not exceed 5. The Lipinski’s parameters were obtained by using the SwissADME server (www.swissadme.ch/index.php) (Daina et al., 2017).

Structure preparation of the proteins and ligands

The crystal structures of the SARS-CoV-2 Spro (PDB ID: 6M0J) and SARS-CoV-2 Mpro (PDB ID: 6Y84) were obtained from the RCSB protein data bank. All the non-standard residues, including water, were removed from the PDB file using Chimera (Pettersen et al., 2004). The RBD domain was obtained from the PDB file of SARS-Cov-2 Spro. The RBD Spro residue sequence number before and after removing ACE2 is presented in . The 3D conformers of the ligands were obtained from the PubChem and optimized using the steepest descent and conjugate gradient steps with General Amber Force Field (GAFF) (Wang et al., 2004) in Chimera (Pettersen et al., 2004).

Figure 1.

Predicted lipophilicity (Log P) values of the spice molecules obtained from different calculation models.

Molecular docking study

The prepared structures of the protein and ligand were subjected to molecular docking analysis using AutoDock Vina (Trott & Olson, 2010). AutoDock Vina is the newest member of the AutoDock family that has improved speed and accuracy. It uses a hybrid scoring function and a quasi-Newtonian optimization algorithm to find the lowest energy confirmations within the search space. A grid box of 40 Å × 65 Å × 70 Å was built with the centre of the box at (11.98, 0.60, 4.79) for the SARS-CoV-2 Mpro. A grid box of size 30 Å × 45 Å × 30 Å with centre at (−36.51, 30.69, 5.48) was prepared for the SARS-CoV-2 RBD Spro. The exhaustiveness of search was set at 20 and 8 for the SARS-CoV-2 Mpro and the SARS-CoV-2 RBD Spro, respectively, to compensate for the larger box volume and reliable results. The docked poses were ranked as per their binding affinities at the end of the docking run. The ligand interactions of the best-docked poses at the active sites of the macromolecule were extracted using PyMol (Schrödinger LLC, 2017). The ligand interactions were analysed using the 2D interaction plot in the Discovery Studio Visualizer (2005). The Coulombic electrostatic potential surface was determined with the help of the APBS plugin available in PyMol (Schrödinger LLC, 2017).

Molecular dynamics (MD) simulation study

To verify the stability of the complex and interaction dynamics of SARS-CoV-2 RBD Spro and SARS-CoV-2 Mpro, we performed the MD simulation study of the two complexes with the highest docking score (Spro-Piperine and Mpro-Piperine), using GROMACS-5.1.5 (Abraham et al., 2015). The CHARMM36 forcefield (Huang & MacKerell, 2013) was used for the simulation of the systems. The topology parameters for the ligand molecule were obtained from CGenFF (Vanommeslaeghe & MacKerell, 2012). A dodecahedron simulation box filled with TIP3P water model (Price & Brooks, 2004) was prepared. Counter ions were added to maintain the electrical neutrality of the system. The systems were kept at a buffer concentration of 0.15 M. Then, the build systems were energy minimized with 50,000 steps using the steepest descent algorithm. Then the systems were equilibrated under NVT and NPT ensembles for 100 ps at 300 K temperature and 1 atm pressure before the production run. After the equilibration, a production run of 100 ns was incorporated under the NPT ensemble. For long-range electrostatic interaction, particle mesh Ewald (PME) (Darden et al., 1993) and for van der Waals interactions, the force-switching scheme was used. Besides, for temperature and pressure coupling, the Berendsen thermostat (Berendsen et al., 1984) with velocity rescaling and Parrinello-Rahman barostat (Parrinello & Rahman, 1981) with isotropic rescaling were used, respectively. The simulation time step was set to 2 fs, and the trajectories were recorded at every 10 ps. The simulation data were analysed by analysing the root mean square deviation (RMSD), root mean square fluctuation (RMSF), number of hydrogen bonds and radius of gyration (Rg) using Gromacs analysis tools. The principal component analysis (PCA) was performed using the g_covar, g_anaeig and g_sham tools of Gromacs. The data were exported to origin 9.0 and plotted for further analysis purposes.

Mm/PBSA binding free energy calculation

The method of calculation of binding free energy from MD trajectory snapshots using the molecular mechanics Poisson–Boltzmann surface area method is widely used. The binding free energy of the systems was estimated by extracting the snaps from the last 20 ns of the MD simulation using g_mmpbsa tool of Gromacs (Baker et al., 2001; Kumari et al., 2014). The binding free energy takes the contribution from vacuum potential energy, polar solvation energy and non-polar solvation energy. The binding free energy can be represented as where Gcomplex, Gprotein and Gligand are the total free energies of the complex, isolated protein and isolated ligand, respectively. The free energy of the individual terms was estimated by where x is the complex, protein or ligand, and TS represents the entropic contribution to free energy in a vacuum with T and S as temperature and entropy. The average molecular mechanics potential and solvation free energies were calculated by using Equations (3) and (4) where Ebonded takes the contribution from a bond, angle and dihedral terms and Enonbonded consists of electrostatic and van der Waals energy contributions. The solvation energy includes the polar and non-polar solvation energies from the Poisson–Boltzmann equation and solvent accessible surface area (SASA), respectively.

Results and discussion

Molecular electrostatic potential (MEP) surface analysis

The electrostatic potential is an essential property for the review and prediction of the reactive behaviour of a molecule. The study of the MEP surface can provide information about the active site of the macromolecule with the indication of relative ligand orientation and nature of the active site at which an approaching electrophile is attracted (Politzer et al., 1985). In a biological macromolecule, the electrostatic potential surface is plotted by analysing the electron-rich and deficient regions of the molecule. The detailed insight at the molecular label helps to predict the potentiality of the ligands to take part in chemical reactions and their mechanism of interaction. The MEP surface representation of SARS-CoV-2 RBD Spro and Mpro ligand-binding sites with the simultaneous presence of all docked molecules are provided in . The electronically poor regions (blue) are referred to as positive potential, whereas the dense electron regions (red) are at a negative potential, and the white zones are considered as neutral. From the figure, it is observed that most of the selected molecules actively bound at the red regions that are referred to as highly negative electrostatic potential regions. This implies that the molecules are polar and can actively take part in the binding process with stable interactions, which in turn could help to block the viral replication.

Figure 2.

Predicted solubility (Log S) values of the spice molecules obtained from different calculation models.

Lipophilicity

Lipophilicity or fat friendliness of a molecule defines the dissolving capability in fat, oil or any non-polar solvent (Lindsley, 2010). The water n-octanol partition coefficient (log Po/w) is used as the measure of lipophilicity (Constantinescu et al., 2019). Various computational methods are developed for the estimation of log Po/w for diverse performance upon different chemical sets. The SwissADME provides five different predictive models such as XLOGP3 (Cheng et al., 2007), WLOGP (Wildman & Crippen, 1999), MLOGP (Moriguchi et al., 1992), SILICOS-IT (http://silicos-it.be.s3-website-eu-west-1.amazonaws.com/software/filter-it/1.0.2/filter-it.html, 2016) and iLOGP (Daina et al., 2014) for better prediction accuracy.

Predicted lipophilicity (Log P) values of the spice molecules obtained from different calculation models are shown in Figure 1. All the molecules subjected to lipophilicity test lie in the range of +1.2 to +4.19 of consensus value that obeys the Lipinski’s limit of log p < 5, which suggests they can be used for further clinical trials (Arnott & Planey, 2012). The lowest lipophilicity is observed for Vanillin and the highest for Nerolidol, among the screened molecules. From Figure 1, it is found that all the ligand molecules have positive lipophilicity value. Hence, these molecules satisfy the essential criteria to be drug molecules.

Water solubility

Solubility (Log S) is the measure of homogeneity of the system from the mixture of solute and solvent. It is considered as one of the vital parameters in drug concentration determination for a desired pharmacological response (Savjani et al., 2012). Poor solubility of drugs is a major issue in drug discovery and development. Solubility acts as a driving force to attain high drug concentration in blood for therapeutic effectiveness (Bergström & Larsson, 2018).

The drug solubility property of the proposed small molecules was obtained from the SwissADME. The server used three solubility models, such as Ali (Ali et al., 2012), ESOL (Delaney, 2004) and Silicos-IT (http://silicos-it.be.s3-website-eu-west-1.amazonaws.com/software/filter-it/1.0.2/filter-it.html, 2016) that is comprised of different topological methods to check the water solubility of these small molecules. The plot for solubility of the proposed small molecules based on these three different models is shown in Figure 2. The Log S values obtained for the ligand molecules based on these three models are in the range of −1.8 to −3.94 for ESOL method, −1.68 to −4.99 for Ali method and −1.48 to −5.52 for Silicos-IT method. The values from different models suggest to the moderately soluble to very soluble nature of the molecules. The reference values of Log S for moderately soluble and highly soluble molecules range from −4 to −6 and −2 to −4, respectively. The solubility values suggest for the oral administration of these molecules.

Pharmacokinetic properties

The pharmacokinetic property is the prime factor for the selection of a drug candidate that describes the drug disposition in the body. The significant parameters that quantify the pharmacokinetics of a drug are its ADME (absorption, distribution, metabolism, excretion) properties (Jang et al., 2001). All the molecules subjected to ADME tests are qualified for drug approval with their high value of gastrointestinal (GI) absorption (Daina & Zoete, 2016), which in turn implies for their use as an oral drug. Table 1 represents the pharmacokinetic properties of the proposed drug candidates. The passive GI absorption and blood–brain barrier (BBB) permeation is a fundamental criterion for the distribution of the drug molecules. From Table 1, it is observed that all the ligand molecules are BBB permeant that implies their underlying distribution index. The high negative skin permeable coefficient (Kp) values indicate a less skin permeability that is useful for their transdermal delivery. The interaction of the drug molecules with cytochromes P450 (CYP) is an essential property as they play a crucial role in drug elimination through biotransformation metabolism. The noninhibition of CYP isoforms such as CYP1A2, CYP2C19, CYP2C9, CYP2D6 and CYP3A4 disclose that these molecules are not the substrate for these enzymes that resembles for the lower degradation rate of these molecules, which will make it effectively available for blocking the SARS-CoV-2 RBD Spro and SARS-CoV-2 Mpro. The synthetic accessibility values suggest the facile synthesis of these molecules. All these parameters infer these close to drug-like molecules, which may be used as successful drug candidates.

Docking study of SARS-CoV-2 receptor-binding domain spike protein

In addition to the above investigations, a molecular docking study was performed to estimate the binding affinity and their binding pose of the ligand molecules at the binding site of the SARS-CoV-2 RBD Spro. From the study, it is observed that Piperine has the highest interaction affinity among the screened compounds. The docked poses of the four ligand molecules (Piperine, Capsaicin, Gingerol and Terpinen-4-ol) along with their 2D interaction diagram having the highest binding affinity, among the selected molecules, are presented in descending order in Figure 3. From Table 1, it is observed that these four molecules follow the trend for their binding affinity with Piperine (−6.4 kcal/mol) at the highest, then Capsaicin, Ginerol and Terpinen-4-ol (all having −5.5 kcal/mol) among all the selected molecules. From Figure 3(a), it is observed that Piperine is associated with hydrogen bond interaction with GLY164 and GLY170. TYR173 (TYR505) and SER162 (SER494) are involved with pi–pi T-shaped and carbon-hydrogen bond interactions, respectively. The binding process is also governed by van der Waals interactions with the residues ARG71, TYR121 (TYR453), TYR163 (TYR495) and ASN169 (ASN501) of SARS-CoV-2 RBD Spro. Hence, the interaction of Piperine with SARS-CoV-2 RBD Spro is stabilized by covalent hydrogen bonding, pi–pi T-shaped and van der Waals interactions with a good affinity score. Capsaicin interacts with the residues GLY164 and TYR173 (TYR505) through pi-Donor hydrogen bond and pi–pi T-shaped interactions with the benzene ring, respectively (Figure 3(b)). The residues ARG71 (ARG403), ASP73, GLU74, GLN77, LYS85, TYR121 (TYR453), SER162 (SER494), TYR163 (TYR495) and ASN169 (ASN501) are involved with van der Waals interaction with Capsaicin.

Lowest energy docked pose of (a) Piperine, (b) Capsaicin, (c) Gingerol and (d) Terpinen-4-ol with SARS-Cov-2 RBD Spro and their 2D interaction diagram. The colour codes represent the nature of interactions.

On the other hand, Gingerol is stabilized by various kinds of interactions with the SARS-CoV-2 RBD Spro (Figure 3(c)). The residues GLY164, ASN169 (ASN501) and GLY170 are associated with hydrogen bond interaction with Gingerol. Other than the hydrogen bond interaction TYR173 (TYR505) is having a pi–pi T-shaped interaction with the benzene ring of Gingerol while ARG71 (ARG403), TYR121 (TYR453), TYR163 (TYR495), PHE165 (PHE497) and GLN166 (GLN498) residues are involved in van der Waals interactions. The Terpinen-4-ol is stabilized by hydrophobic interaction with the residues ARG125, LYS126, TYR141 and PRO159 while the residues ARG122, PHE124, ASP135, SER137, GLU139 and ILE140 are involved in van der Waals interactions with Terpinen-4-ol (Figure 3(d)). The lowest energy poses of the rest 26 molecules along with their 2D interaction diagrams are provided in .

Docking study of SARS-CoV-2 main protease

The above-selected molecules were also docked with the SARS-CoV-2 Mpro to observe the inhibitory effect of these molecules. The docking study reveals that all the molecules are interacting with the SARS-CoV-2 Mpro with certain binding affinity. The docking data are also presented in Table 1. From Table 1, it is seen that Piperine has the highest affinity at the binding site of SARS-CoV-2 Mpro among all the selected molecules, which is similar to the case of SARS-CoV-2 RBD Spro. The ΔG value, known as binding free energy, for the four molecules having the highest affinity among all the selected molecules, along with their 2D interaction diagram, is given in Figure 4. The four molecules have followed the binding affinity trend as Piperine (−7.3 kcal/mol) > Capsaicin (−6.4 kcal/mol) > Carvone (−6.2 kcal/mol) > Gingerol (−6.1 kcal/mol). From Figure 4(a), it is observed that the interaction of Piperine at the binding site of the SARS-CoV-2 Mpro is stabilized by hydrogen bonding, electrostatics and van der Waals interactions. The residues GLN299 and VAL303 are associated with hydrogen bonding interaction; ASP295 and ARG298 with charged interactions while MET6 and PRO9 are associated with hydrophobic interactions with Piperine. The molecule is also stabilized through van der Waals interactions with residues PHE8, GLY127, ILE152, PHE291 and THR304 at the binding site of the SARS-CoV-2 Mpro. The molecule Capsaicin is stabilized in the binding pocket through van der Waals and hydrophobic interactions (Figure 4(b)). The residues MET6, PHE8, PRO9 and ILE152 are interacting through hydrophobic interactions such as alkyl and pi–alkyl with the Capsaicin. Capsaicin is interacting with residues ALA7, GLY11, LYS12, GLN127, TYR154, PHE291, ASP295, ARG298, GLN299, VAL303 and THR304 through van der Waals interaction. The interaction of Carvone with the SARS-CoV-2 Mpro is stabilized through hydrophobic and van der Walls interactions (Figure 4(c)). Carvone interacts with the residues MET6, PHE8 and ARG298 of SARS-CoV-2 Mpro through hydrophobic contacts. The residues ALA7, PRO9, GLN127, ASP295, GLN299, GLY302 and VAL303 are in van der Waals interactions with Carvone. Gingerol is stabilized by hydrogen bond, hydrophobic and van der Waals interactions in the binding pocket of the SARS-Cov-2 Mpro (Figure 4(d)). VAL303 is interacting through hydrogen bond interaction with Gingerol. The residues LYS12 and THR304 are involved in carbon-hydrogen bond interactions while MET6, ALA7, PHE8, GLN127, TYR154, ASP295 and ARG298 are associated with van der Waals interactions with Gingerol. The interaction of Gingerol with the residues PRO9 and ILE152 is stabilized through pi–alkyl interactions. The lowest energy binding poses of the rest 26 molecules along with the ligand interaction diagram at the binding sites of the SARS-CoV-2 Mpro are provided in .

Lowest energy docked pose of (a) Piperine, (b) Capsaicin, (c) Carvone and (d) Gingerol with SARS-CoV-2 Mpro and their 2D interaction diagram. The colour codes represent the nature of interactions.

The MD simulation is one of the proven in-silico methods for the determination of protein–ligand dynamics concerning a high temporal resolution of nanosecond or picosecond order. Here the docked poses of RBD Spro and Mpro with Piperine were used for a 100 ns MD simulation to analyse the stability of these docked compounds.

Root mean square deviation (RMSD)

The RMSD values from MD simulation provide information about structural and conformational stability. Figure 5 represents the backbone RMSD data of viral proteins and their complex with Piperine. From the plot, it is observed that both the simulations have less fluctuation throughout the simulation time. The average RMSD values of RBD Spro, RBD Spro-Piperine, Mpro and Mpro-Piperine are calculated as 0.143 ± 0.025 nm, 0.130 ± 0.018 nm, 0.212 ± 0.041 nm and 0.203 ± 0.028 nm, respectively. The average RMSD values of the Piperine bound proteins as compared to only proteins are found to be less representing to their conformational stability. Both the simulations are attained equilibrium within 0.3 nm, which is also a measure of the systems’ stability during the simulation (Al-Shabib et al., 2018, 2020; Millan et al., 2018).

Figure 5.

Root mean square deviation plots of (a) RBD Spro (black) and RBD Spro-Piperine (red); (b) Mpro (black) and Mpro-Piperine (blue).

Root mean square fluctuation (RMSF)

The conformational fluctuations of the proteins RBD Spro and Mpro were analysed by observing the residual changes that resulted due to the interaction of Piperine with the proteins. The RMSF plots of the Cα atoms of the viral proteins and their complex with Piperine are shown in Figure 6. From the analysis, it is found that RBD Spro, RBD Spro-Piperine, Mpro and Mpro-Piperine have the average RMSF values 0.099 ± 0.060 nm, 0.097 ± 0.051 nm, 0.119 ± 0.077 nm and 0.120 ± 0.077 nm, respectively. It is observed that RBD Spro-Piperine (Figure 6(a)) and Mpro-Piperine (Figure 6(b)) show similar fluctuations as compared to only RBD Spro and Mpro, which implies to the stability of these compounds. In addition to that, a majority of the protein residues are found to be stabilized within RMSF 0.3 nm. The decrease in fluctuations of Piperine bound to RBD Spro also suggests for the active binding of Piperine (A. Kumar, Choudhir, et al., 2020).

Figure 6.

Root mean square fluctuation plots of Cα-atoms of (a) RBD Spro (black) and RBD Spro-Piperine (red); (b) Mpro (black) and Mpro-Piperine (blue).

Radius of gyration (Rg)

The root mean square distance between an object and the centre of gravity is defined as the radius of gyration (Rg). The radius of gyration is a measure of the compactness of the protein structure, where higher Rg value is referred to as a less compact structure, and low Rg value is inferred as high compactness that implies more stability. The measured average Rg values of RBD Spro, RBD Spro-Piperine, Mpro and Mpro-Piperine are 1.829 ± 0.008 nm, 1.833 ± 0.010 nm, 2.233 ± 0.012 nm and 2.238 ± 0.013 nm, respectively. From Figure 7(a,b), it is observed that there is a little enhancement in the Rg values of RBD Spro-Piperin, and Mpro-Piperine as compared to RBD Spro and Mpro, which implies to the gain in compactness of the protein structures upon binding to Piperine.

Figure 7.

Radius of gyration plots of (a) RBD Spro (black) and RBD Spro-Piperine (red); (b) Mpro (black) and Mpro-Piperine (blue); (c) Intermolecular hydrogen bonds formed between RBD Spro-Piperine and Mpro-Piperine during 100 ns MD simulation.

Number of hydrogen bonds

The number of hydrogen bonds formed between the protein–ligand complex is the measure of the binding strength of the ligand to the protein. The RBD Spro (red) and Mpro (blue) bound to Piperine have a constant number of 1–2 hydrogen bonds throughout the simulation time (Figure 7(c)). There is a maximum number of 3 and 4 hydrogen bonds observed in the case of RBD Spro-Piperine and Mpro-Piperine, respectively. The number of hydrogen bonds fluctuates throughout the simulation time for both RBD Spro-Piperine and Mpro-Piperine, which suggests for conformational changes in the binding site of the ligand during the simulation. The observation from hydrogen bond analysis indicates that the complexes are stable for the performed simulation time.

Interaction energy

The interaction energy is the measure of the interaction strength of the protein–ligand complex. In order to validate the results of molecular docking studies, the analysis of the interaction free energies from MD simulation was performed. The average interaction energy takes the contribution from the average short-range Lennard-Jones and van der Waals energy. The average interaction energies of RBD Spro-Piperine and Mpro-Piperine are found to be − 41.401 ± 17.843 kJ/mol and −143.162 ± 23.043 kJ/mol, respectively. These interaction energy values suggest that Piperine binds to the RBD Spro and Mpro with good affinity and hence supports the docking results, which in turn helps for the favourable use of Piperine as a drug candidate for SARS-CoV-2.

Solvent accessible surface area (SASA)

SASA is a measure of the receptor exposure to the solvent environment during the simulation. The hydrophobic residues that got exposed to the solvent environment upon binding with the ligand molecules contribute to the SASA values. The plot of the SASA for the proteins and their ligand-bound form is presented in Figure 8. The analysed average SASA values for the RBD Spro, RBD Spro-Piperine, Mpro and Mpro-Piperine are 106.976 ± 1.602 nm2, 107.235 ± 1.667 nm2, 150.698 ± 2.565 nm2 and 151.022 ± 2.207 nm2, respectively. There is no significant change observed for the averaged SASA values of the complex as compared to only protein suggesting their stability after binding to the drug molecule.

Figure 8.

Solvent accessible surface area (SASA) plots of (a) RBD Spro (black) and RBD Spro-Piperine (red); (b) Mpro (black) and Mpro-Piperine (blue).

MMPBSA binding free energy analysis

MD simulation can also be used to calculate the binding free energy of the protein–ligand complex. The binding free energy is the measure of the stability of the system in turns of consistency of nonbonded interactions throughout the simulation. The binding free energy was calculated by using MMPBSA method by taking 2000 snapshots from the trajectory. The computed value of binding free energy for RBD Spro-Piperine is found to be −5.533 ± 0.839 kJ/mol, and for Mpro-Piperine is −37.971 ± 0.271 kJ/mol. It is observed that for both RBD Spro and Mpro, van der Waals energy plays a crucial role in the interaction process. The van der Waals energy, electrostatic energy and non-polar energy are contributed actively to the total interaction energy. In contrast, polar energy has a positive contribution to the whole interaction process. The observed data indicate that the van der Waals, electrostatic and non-polar interactions combinedly contribute to the stability of both the compounds. The contribution from different interactions to the binding free energy for RBD Spro-Piperine and Mpro-Piperine is provided in .

Table 2.

Lowest energy binding affinity of Piperine and few of the currently used drugs for SARS-CoV-2 as obtained from molecular docking study.

Molecule	Binding affinity (kcal/mol)
	RBD Spro	Mpro
Piperine	–6.4	–7.3
Chloroquine	–5.0	–4.9
Favipiravir	–5.3	–5.6
Hydroxychloroquine	–4.8	–6.0
Oseltamivir	–5.1	–5.5
Remdesivir	–6.1	–7.2
Ribavirin	–5.6	–6.1

Principal component analysis (PCA)

The PCA is an essential technique to monitor the conformational dynamics of biomolecules. It is useful in determining the concerted motion of protein as well as protein–ligand complex from the MD trajectories. The diagonalization of the covariance matrix of backbone atoms of the proteins and ligand-bound form were considered for the principal components PC1 and PC2 (Figure 9). From Figure 9(a,b), it is observed that both the Spro-Piperine and Mpro-Piperine are less flexible as compared to unbound proteins since they covered less conformational space. It concludes that the ligand-bound forms are more stable as compared to the unbound proteins.

Principal component analysis of (a) RBD Spro (black), RBD Spro-Piperine (red) and (b) Mpro (black) and Mpro-Piperine (blue). Free energy landscape plot of (c) RBD Spro, (d) RBD Spro-Piperine and (e) Mpro and (f) Mpro-Piperine.

The principal components obtained were used as the reaction coordinates to find the Gibbs free energy landscape (Figure 9) to visualize the energy minima of the unbound protein as well as the protein–ligand complex. From Figure 9(c–f), it is observed that both the ligand-bound proteins have less Gibbs-free energy values than the unbound proteins indicating their stability and energetically favourable conformational transitions. The shape and size of the minimum energy area (blue colour) in case of RBD Spro-Piperine and Mpro-Piperine are more as compared to the unbound proteins RBD Spro and Mpro, which suggests the ligand-bound forms are thermodynamically more favourable.

The comprehensive study reveals that Piperine forms a stable complex with RBD Spro and Mpro and can be considered as an active inhibitor against SARS-CoV-2. From the docking results, it is observed that the Piperine molecule is the best candidate for the inhibition of the RBD Spro and the Mpro of SARS-CoV-2 among the selected 30 molecules. To observe the effectiveness of Piperine over currently used drugs, we carried out the docking study of a few drug molecules such as chloroquine, favipiravir, hydroxychloroquine, oseltamivir, remdesivir and ribavirin using the same docking protocol as followed for the 30 spice molecules. From the docking score, it is found that Piperine performed better as compared to the currently used drugs stated above. The lowest energy pose of a few presently used drugs with their 2D interaction diagram is provided in Supplementary Figures S5 and S6 corresponding to SARS-CoV-2 Mpro and RBD Spro, respectively. A comparison of the lowest energy dock scores of these drug molecules along with Piperine is also provided in Table 2. The MD simulation results reveal that Piperine actively inhibits both the RBD Spro and Mpro by binding to their active sites. Piperine binds on the active site of the RBD Spro with those residues by which it interacts with ACE2. So, the binding of Piperine on that site may potentially cease the interaction tendency of RBD Spro with ACE2. Similarly, the interaction of Piperine on the active site of the Mpro may inhibit its viral replication. From the docking and MD results, we conclude that Piperine forms a very stable complex with RBD Spro and Mpro and shows better affinity as compared to the currently used drugs that are mentioned above against SARS-CoV-2.

Conclusion

This study used molecular docking and MD simulation as potential tools to monitor the inhibitory efficiency of natural spice molecules against SARS-CoV-2, which emerged as a global threat to millions of people across the globe. It is observed that all the proposed spice molecules qualified the ADME test with their suitable pharmacokinetic properties to be useful as a drug candidate. The docking study is revealed that all the molecules actively take part in binding to the SARS-CoV-2 RBD Spro and Mpro with their low or high value of binding affinity. This binding of these molecules will help to inhibit the replication of the viral proteins with specific hindrances upon their mutarotation. For both the viral targets, Piperine performed well with its highest binding affinity of −6.4 and −7.3 kcal/mol for SARS-CoV-2 RBD Spro and Mpro, respectively. Besides, Piperine is also found more effective as compared to a few of the currently used drugs. The MD simulation study is supported for the stable interaction of Piperine with SARS-Cov-2 RBD Spro and Mpro. The PCA and binding free energy results also suggest for the active participation of Piperine in stable complex formation with RBD Spro and MPro. Hence, the study proposes Piperine as an active molecule for the inhibition of SARS-CoV-2. Since this study is performed computationally, it requires wet-lab experiments in-vitroas well as in-vivofor further validation.

Click here for additional data file.