Identification of new anti-nCoV drug chemical compounds from Indian spices exploiting SARS-CoV-2 main protease as target.

The 2019-novel coronavirus (nCoV) has caused a global health crisis by causing coronavirus disease-19 (COVID-19) pandemic in the human population. The unavailability of specific vaccines and anti-viral drug for nCoV, science demands sincere efforts in the field of drug design and discovery for COVID-19. The novel coronavirus main protease (SARS-CoV-2 Mpro) play a crucial role during the disease propagation, and hence SARS-CoV-2 Mpro represents as a drug target for the drug discovery. Herein, we have applied bioinformatics approach for screening of chemical compounds from Indian spices as potent inhibitors of SARS-CoV-2 main protease (PDBID: 6Y84). The structure files of Indian spices chemical compounds were taken from PubChem database or Zinc database and screened by molecular docking, by using AutoDock-4.2, MGLTools-1.5.6, Raccoon virtual screening tools. Top 04 hits based on their highest binding affinity were analyzed. Carnosol exhibited highest binding affinity -8.2 Kcal/mol and strong and stable interactions with the amino acid residues present on the active site of SARS-CoV-2 Mpro. Arjunglucoside-I (-7.88 Kcal/mol) and Rosmanol (-7.99 Kcal/mol) also showed a strong and stable binding affinity with favourable ADME properties. These compounds on MD simulations for 50 ns shows strong hydrogen-bonding interactions with the protein active site and remains stable inside the active site. Our virtual screening results suggest that these small chemical molecules can be used as potential inhibitors against SARS-CoV-2 Mpro and may have an anti-viral effect on nCoV. However, further validation and investigation of these inhibitors against SARS-CoV-2 main protease are needed to claim their candidacy for clinical trials.Communicated by Ramaswamy H. Sarma.

Introduction

Since December 2019, a outbreak of COVID-19 with massive global impact has started in Hubei and Wuhan city in China caused by a novel coronavirus, SARS-CoV-2 (Bhoopathi et al., 2020; Kumar & Rathi, 2020; Wang et al., 2020). On January 2020, World health organization emergency committee declared a global health emergency based on the high rate of spreading of the infection with high fatality rate (Chhikara et al., 2020; Kumar & Rathi, 2020). As of now on 4th April 2020, Europe and USA are new epicentres of COVID-19 Pandemic. Even in the past, different coronaviruses have caused multiple human diseases that resulted in global epidemics such as the Middle East respiratory syndrome (MERS), severe acute respiratory syndrome (SARS) and coronavirus disease 2019 (COVID-19). While the coronaviruses have a significant impact on human health, the general public has an inferior awareness of coronavirus pathogenesis and infection. COVID-19 was declared a pandemic in March 2020 as the worldwide human population is facing high risk from contracting the nCoV infection. Many countries, including India, has announced lockdown in the country and to maintain social distancing to avoid further spread as no drug or vaccine is available against SARS-CoV-2. Coronavirus has been classified into four subfamilies based on their shape and host. These subfamilies are alpha-coronavirus, beta-coronavirus, gamma-coronavirus and delta coronavirus (Paules et al., 2020). Alpha- coronavirus and beta-coronavirus are considered to have been originated from bats, while the gamma and delta coronaviruses are considered to be derived from birds and pigs (Banerjee et al., 2019). Coronaviruses contain a positive sense, single-strand RNA genome coding for viral polymerase, RNA synthesis materials, and large non-structural polypeptide. Coronavirus genome contains transcriptional modification, including 5’methylated cap and a 3’polyadenylated tail. Coronaviruses have very high rates of error in RNA replication due to constant errors by RNA dependent RNA polymerase (Banerjee et al., 2019; Lau et al., 2018).

Only a few protein crystal structures of SARS-CoV-2 are available on protein databank. SARS-CoV-2 main protease, a potential drug target, crystal structure (PDB-ID: 6Y84) was available and used for docking simulation and identification of potential drug molecule form Indian spices. SARS-CoV-2 main protease has a vital role in the processing of polyprotein that is translated from viral RNA, and the protease is considered as key for viral survival and growth (Prajapat et al., 2020).

Medicinal plants yielding biologically active compounds have always been of great interest to scientists as they play an essential role in preventing human diseases. In the entire world, India is recognized where spices have been traditionally used as a source of medicine. Many active pharmaceutical gradients have been identified and extracted having a wide range of physiological and pharmacological properties (Sachan et al., 2018). Apart from the regular uses of spices in culinary activities, they are widely used as indigenous medicines, nutraceuticals, in aromatherapy, as natural colouring agents, perfumes, cosmetics etc. In recent years there has been experimental evidence on physiological benefits that could be drawn in the context of various diseases like diabetes, cardiovascular issues and inflammatory disorders like arthritis and cancer. They have also been identified as preventive agents in certain conditions. Spices like red pepper, garlic, and fenugreek have been reported to have hypercholesterolemic activity, whereas fenugreek and garlic are also known to reduce and control blood sugar levels (Bhagya & Raveendra, 2017). The main active ingredient of turmeric, Curcumin has been widely studied to have a broad spectrum of medicinal value ranging from anti-cancerous, anti-inflammatory and an anti-amyloidogenic activity (Bhagya & Raveendra, 2017). Several computational studies related to drug or vaccine development against SARS-CoV-2 is recently published (Aanouz et al., 2020; Gupta et al., 2020; Joshi et al., 2020; Muralidharan et al., 2020; Sarma et al., 2020).

In the current study, we have utilized the prior knowledge on the medicinal values and potential applications of Indian spices, we have tried to explore if they can be used as novel agents for controlling SARS-CoV-2. The data generated using computational approaches give very encouraging results.

Material and methods

Data sources and preparation target protein

In this study, prepared list of 45 chemical compounds from Indian spices and compound structure file was downloaded from PubChem database or Zinc database. Structure of protein SARS-CoV-2 Mpro was downloaded from RCSB protein database. (PDB ID: 6Y84).

The atomic coordinates of active site were defined using report available in the literature and identification of catalytic His41, His164 and Cys145 (Khaerunnisa et al., 2020). Energy minimization of the protein molecule was done by Swiss PDB Viewer (Guex & Peitsch, 1997). Before docking, the assignment of charge, solvation parameters and fragmental volumes to the protein was done using the Autodock Tool 4 (ADT). The protein 6Y84 PDB molecule was further optimized by using ADT for the molecular docking using established procedure.

The structure data files of the all chemical compounds were downloaded from PubChem database and converted into mol2 structures by using Open Babel. In order to further simplify the analyses, ligands were fist optimized and converted mol2 to PDBQT format by using the graphical user interface version of Raccoon. MGLTools-1.5.6, Raccoon is preparing AutoDock virtual screening tool-python (Forli, Scripps Research Institute).

Compound screening using raccoon |AutoDock program

Molecular screening of the all chemical compound libraries was performed by using Raccoon and MGLTools-1.5.6 software by Autodock as the engine for docking (Morris et al., 2009). During the molecular docking period, the ligands were considered as a flexible molecule and the protein was considered as a rigid structure molecule. The configuration file for the grid parameters file and docking parameters file was generated by using Autodock. Autodock and Autogrid tools integrated with the Autodock4 were used to generate grid maps for each atom of the ligand. The grid boxes were made, such as to include one site at a time and perform docking. Grid X, Y and Z coordinates were 9.204, −4.557, and 19.602. We analyzed each ligand by setting default docking parameters except in the number of runs: We ran the LGA for 100 runs with each ligand with an initial population size of 150 random starting positions and conformation, 2.5 million number of energy evaluations.

The application of grid parameters file was also used to predict the amino acid residue in the active site of the protein that interacts with the ligands. Positional root-mean-square deviation (RMSD) values were less than 1.0 Å considered ideal and clustered together for finding the favourable binding. The highest negative binding energy was considered as the ligand molecule with maximum binding affinity.

Analyses and visualization

Visual analysis of the docking site was performed using Pymol-2.3.3 and the results were validated by using AutoDock Tools-1.5.6 (Morris et al., 1998). Binding interaction analyses of identifying inhibitor and the SARS-CoV-2 Mpro (PDB ID: 6Y84) was done using an online program by using Ligplot analysis (Wallace et al., 1995).

MD simulations

The classical molecular dynamics simulation is carried out for the prepared co-crystal structure and selected ligand-binding complex poses through Desmond Molecular dynamics package incorporated in Schrodinger suite (Chow et al., 2008). Both the apo and ligand complex is solvated in the TIP3P model (specifies a three-site rigid water molecule with charges), using the volume occupancy in an orthorhombic box with periodic boundary conditions. For neutralizing the system, the overall charge of apo and ligand complex is solvated with appropriate cation (Na+) or anion (Cl-) along with a salt concentration of 0.15 mol/L. The prepared system is energy minimized for a convergence threshold of 1.0 kcal/mol/Å by using the steepest descent method, and for minimization and relaxation of the system, the NPT ensemble is applied (Selvaraj et al., 2015). The standard temperature is kept constant at 300 K and pressure at the level of 1.013 bar for the total simulation, and each simulation is started for the total time scale of 50 ns. For the analysis of MD trajectories, the RMSD, RMSF and Hydrogen bond interactions are analyzed using the trajectory analysis incorporated in the Desmond (Selvaraj & Singh, 2014).

ADME analyses

The absorption, distribution, metabolism, and excretion (ADME) properties of the studied top hits compound were calculated by using online SwissADME Program (Guex & Peitsch, 1997). The significant parameters for ADME associated properties such as Lipinski’s rule of five, pharmacokinetic properties the solubilities of drug and drug likeness were considered.

Result and discussion

We have created a database of 45 Indian spices compounds and their structures downloaded from PubChem and Zinc databases. Four small molecules were selected based on their binding affinity with SARS-CoV-2 Mpro as shown in the Table 1. The list of remaining compounds used of moleculer docking is submitted as supplementary material (Table 1S). The top three compounds include Carnosol, Rosmanol, and Arjunglucoside-I. The further insights into binding interactions of these compounds with SARS-CoV-2 Mpro were analysed using ligplot as shown in Figures 1 and 2.

Table 1.

Details of various kinds of interaction shown between the amino acids near the active site of SARS- CoV-2 main protease along with their respective inhibitor constant (Ki) and biological source and binding energy.

Compound	Source of the compound	Binding affinity (kcal/mol)	Amino acid residues	Inhibitor Constant (Ki)
Alpha-ketoamide	Positive Control	−9.48		112.59 nM
Carnosol	R. officinalis	−8.2	Cys145, His164, Glu166, Gln189, Met165, Arg188, Cys44, Met49, *His41, Thr45, Thr25, Thr26, Leu27, Gly143, Ser46, Asn142	969.58 nM
Rosmanol	R. Officinalis	−7.99	Cys145, Gly143, Thr25, His41, Met49, Asn142, His143, Glu166, Leu141, Ser144, Met165, Phe140	1.38 µM
Arjunglucoside-I	Terminalia chebula	−7.88	Cys145, Ser144, Leu141, Met165, Phe140, His163, Glu166, Asn142, His41, Gly143, Thr26, Leu27, Thr25	1.67 µM

The active site residues are indicated in bold.

Figure 1.

Ligplot images showing both hydrogen and hydrophobic interactions by (a) α Ketoamide and (b) Carnosol with SARS- CoV-2 main protease. Circled residues represent interaction with active site residues. Cysteine145 forms hydrogen bonding with both the compounds.

Figure 2.

Ligplot images showing both hydrogen and hydrophobic interactions by (a) Rosmanol and (b) Arjunglucoside with SARS- CoV-2 main protease. Circled residues represent interaction with active site residues.

SARS-CoV-2 Mpro Protein (PDBID- 6Y84) active site defined but unliganded (Owen et al., 2019) was used and the molecules where docked in the region where α–ketoamide was bound to protein (PDBID- 6Y2F) (SARS- CoV-2) main protease with bound with α–ketoamide (Zhang et al., 2020) involving His41,164 and Cys145 in the active site. α–ketoamide was hence used as a positive control for this study and resulting interactions were analytically compared with the other Indian spices. α–ketoamide showed interactions with Cys145 via hydrogen bonds and hydrophobic interactions His41. It was also seen to form hydrogen bonds with Thr26. Upon further analysis of the docking results, it was seen that our molecules had comparable binding energy as to α–ketoamide suggesting that these ingredients could indeed interact with same site (amino acids) as α–ketoamide. Carnasol (-8.2 Kcal/mol) was seen to form hydrogen bonds with Leu141, Ser144 and Cys145. Hydrophobic interactions for Carnasol include His41, Thr25, Asn142, Phe140, Glu166, Met165 (Figure 1). Rosmanol, isolated from Rosemary (Rosmarinus officinalis), has been previously reported to have antioxidant activity, shows binding energy of -7.99 Kcal/mol and forms hydrogen bonds with Leu141, Gly 143, Ser144 and Cys145. Further hydrophobic interactions are shown with Thr25, 26, His41, Phe140, His163, and Leu27 (Figure 2a). Arjunglucoside-I binds with hydrophobic interaction with His41, 164 and Cys145 additional to hydrogen bonding with Thr25, 26, His163, and Glu166 (Figure 2b).

MD simulations analysis

The results of MD simulation for both apo and ligand complex is analyzed for the 50 ns of time scale to understand the dynamic behavior and stability. MD simulation is performed for the total of 50 ns and the trajectories are for RMSD plot as shown in the Figure 3. The figure represents the colors with (a) Apo protein, complexed with (b) Alpha ketoamide, (c) Arjun glycoside, (d) Carnosol, and (e) Rosmanol for the 50 ns of time scale. For the apo protein, observed the initial 45 ns of MD simulation as stable, but after the 45th ns there is a sudden drift and that may occur due to the higher occupancy of flexible loops in the C-terminal region (182-304 residues). Even though the apo protein shows the value of ∼2.0 Å till 47th ns and remaining 3 ns ends with ∼3.4 Å. For understanding the residue wise fluctuation between the apo and ligand complex, the RMSF values are plotted. The RMSD deviation for the apo and complex proteins clearly shows the dynamic movement that occurs in the loop regions. For understanding the residues participated in the causative for fluctuations, the RMSF plot is provided in the Figure 4. Here, the amino acid in each position is calculated for its deviation value with the calculated timescale of 50 ns. For apo protein, the amino acid from 70-75th, 160-175th position shows the deviation up to ∼4Å, apart from these regions other amino acids are deviating from ∼1-3Å. The deviation that occurs in the 70-75th, 160-175th position may be functionally reason for the drift in apo protein RMSD at 45th ns. While comparing the apo RMSF values with the complex protein RMSF, notably Carnosol is showing high deviations than the other ligand complexes. The amino acid positions from 40-60th, 150-155th, 178-200th position is showing higher deviations, that ranges from ∼4-6Å deviations. Due to these positional amino acid fluctuations, the RMSD of Carnosol deviates approximately 4 Å after the 26th ns. Similarly, another compound Rosmanol also showing the fluctuations in C-terminal regions, that results in the ∼3-3.5 Å deviations in the RMSD values. The protein secondary structure is framed as 3 sheets, 7 beta hairpins, 9 beta bulges, 13 strands, 32 beta turns, 3 gamma turns from this 182-304 residues are occupied dominantly by loop regions. These 122 residues in the C-terminal functionally act in the MD simulation, and thus the sudden drift happens in the 45th ns. This values of RMSD for the apo protein, is compared with ligand molecules in the Figure 3, and for understanding the deviations, the 25th to 50th ns is focused. The results of ligand complex for the MD simulation of 50 ns of timescale shows that the ligand complex Alpha ketoamide (Red color) and Arjun glycoside (Blue color) shows stable throughout the simulation. While comparing these two Alpha ketoamide and Arjun glycoside ligand complex with apo protein, the ligand complex is matched with apo protein for the timescale of 45th ns. After the 45th ns the apo protein is drifted upwards, but the Alpha ketoamide and Arjun glycoside ligand complex remains stable. This may be due to the strong interaction pattern seen in Alpha ketoamide and Arjun glycoside interactions with the protein. Both the Alpha ketoamide and Arjun glycoside ligand bound complex are stable and positioned in the range of ∼2Å, which is close to the apo protein till 47th ns. The Figures 5a and 5b shows the hydrogen bond interactions for the Alpha ketoamide and Arjun glycoside, which clearly shows the minimum participation of 2-3 hydrogen bonds seen in between the Alpha ketoamide and protein, and minimum participation of 4-6 hydrogen bonds between the Arjun glycoside and protein. This active participation of hydrogen bonds between the Alpha ketoamide and Arjun glycoside with protein makes the ligand complex stable for 50 ns of MD simulations. The ligand Carnosol shows stable movement in terms of RMSD values by showing a narrow graph, but while comparing with apo protein, the Carnosol ligand complex is deviated from the 14th ns and stays in the range bound of ∼3.4 Å till the end of the simulations. For attaining the stable MD simulation for the Carnosol ligand bound complex, 2-3 hydrogen bonds are actively contributed and makes the complex stable in the dynamic state. The ligand Rosmanol bound complex shows stability till the 26th ns with the range of ∼2.5 Å, but after that, the drift happens to make the RMSD value deviated in upward direction and fix the positional RMSD with the range of ∼4 Å. This may be due to the loss of hydrogen bonding interactions after the 27th ns seen in the Figure 5d. Figure 5d says that the initial 27 ns shows the hydrogen bond interactions in the range of 2-3 between the protein and ligand, but after 27th ns the ligand losses the hydrogen bonding ability and shows the hydrogen bond interactions in the range of 1-2 between the protein and ligand. Overall, the apo protein shows a narrow range of stability and notable fluctuations are also shown, that indicates the participation of loop structures. Those fluctuations are arrested through the active interactions of ligand molecules that shows strong binding between the protein and ligand.

Figure 3.

RMSD analysis of Apo and ligand complex (C-Alpha) in molecular dynamics simulations for the time scale of 50 ns. The colors represented in the figure (a) Apo protein, (b) Alpha ketoamide, (c) Arjunglucoside, (d) Carnosol, and (e) Rosmanol.

Figure 4.

RMSF analysis of Apo and ligand complex in molecular dynamics simulations for the time scale of 50 ns shows the deviation of each amino acid positions.

Figure 5.

Hydrogen bonds interaction between the protein and the ligand molecules. (a) Alpha ketoamide, (b) Arjunglucoside, (c) Carnosol, and (d) Rosmanol.

We have also analysed the molecules we report here for violation of Lipinski’s rule (Table 2). Rosmanol, Carnosol fits perfectly as within the defined parameters for non-violation of Lipinski’s rule. The molecules have Log P values ranging from 1.15 to 3.27 which imply that these can effectively have suitable cell membrane permeability.

Table 2.

Lipinski’s parameters for drug likeliness and ADMET properties of chosen ligands along with the standard drugs approved.

Molecule	LogP	TPSA (Angstrom)	MW (Da) (g/mol)	H-Donor	H Acceptor	Log S	Violations
Arjunglucoside-I	1.15	197.37	666.84	8	11	−2.09	3
Carnosol	2.97	66.76	330.42	2	4	−4.77	0
Rosmanol	2.50	86.99	346.42	3	5	−3.64	0

Log P value represents the lipophilicity of the molecule whereas Log S value represents the water solubility. TPSA is Total Polar Solvent Accessibility. BA is oral bioavailability of drug and BBB is blood brain barrier permeability. Violation value land between 0-5, in Lipinski’s rule of five, if value is 0, then Drug likeness shows yes and if violations >0 then shows No.

Their number of hydrogen bond donors as well acceptors are well within range for Carnosol and Rosmanol but Arjunglucoside-I show high number of hydrogen bond donors and acceptors. Arjunglucoside is a large molecule having high molecular weight and total polar solvent area. Despite these factors it might prove to be essential in terms of potential drug once preceded with advanced studies.

Conclusion

Amidst the unforeseeable outbreak of CoVID-19 there has been a sudden rise in demand of drug development, vaccines and identification of potential bioactive molecules which could prove to be useful fulfilling the purpose of broadening treatment options. In the quest for finding novel treatment regimen for these kinds of viral outbreaks, screening of already known molecules could also prove to be vital. In this context, we report here some active pharmaceutical ingredients which are present in the commonly used spices in India could prove to be useful. Preliminary in silico investigations show that indeed some molecules like Carnosol and Rosmanol have the properties which can further exploited and investigated for drug candidate against SARS-CoV-2. Although it is imperative to understand that development of rigid and highly specific treatment options will require further expeimental studies.