Systematic virtual screening in search of SARS CoV-2 inhibitors against spike glycoprotein: pharmacophore screening, molecular docking, ADMET analysis and MD simulations.

In the absence of efficient anti-viral medications, the coronavirus disease 2019 (COVID-19), stemming from severe acute respiratory syndrome coronavirus-2 (SARS CoV-2), has spawned a worldwide catastrophe and global emergency. Amidst several anti-viral targets of COVID-19, spike glycoprotein has been recognized as an essential target for the viral entry into the host cell. In the search of effective SARS CoV-2 inhibitors acting against spike glycoprotein, the virtual screening of 175,851 ligands from the 2020.1 Asinex BioDesign library has been performed using in silico tools like SiteMap analysis, pharmacophore-based screening, molecular docking using different levels of precision, such as high throughput virtual screening, standard precision and extra precision, followed by absorption, distribution, metabolism, excretion and toxicity analysis, and molecular dynamics (MD) simulation. Following a molecular docking study, seventeen molecules (with a docking score of less than - 6.0) were identified having the substantial interactions with the catalytic amino acid and nucleic acid residues of spike glycoprotein at the binding site. In investigations using MD simulations for 10 ns, the hit molecules (1 and 2) showed adequate compactness and uniqueness, as well as satisfactory stability. These computational research findings have offered a key starting point in the field of design and development of novel SARS CoV-2 entry inhibitors with appropriate drug likeliness.

Introduction

The novel coronavirus disease 2019 (COVID-19), caused by deadly, pathogenic and infectious novel coronaviruses (nCoV-19), has caused an unprecedented pandemic worldwide due to increased virulence, morbidity and rapid spread beyond the previous severe acute respiratory syndrome (SARS) pandemic of 2002–03 [1, 2]. The very first outbreak was reported on December 31, 2019, and suspected to be assembled in the laboratory of Wuhan, China [3-6]. Citing the exceedingly concerning situation, COVID-19 has been designated as an emergency and a worldwide pandemic by the World Health Organization (WHO) on March 11, 2020 [7]. According to the latest WHO estimates, there were 326,279,424 confirmed cases of COVID-19 worldwide as of January 17, 2022, with 5,536,609 fatalities [8]. The severe effects of quarantine-related counter-measures imposed by ruling governments, such as lock-down, sanitization of public premises, social distancing, closure of academies or sports clubs, cancellation of public or social events, etc., have severely hampered normal life around the world. As a result, the socioeconomic consequences of the epidemic have steadily damaged societal health, education, economy and several vital sectors have been greatly challenged or ceased [9, 10].

Coronaviruses (CoVs) are single-strand positive-sense RNA enclosed viruses that belong to a large family [11]. Being club-shaped glycoprotein, they are around 60–140 nm in diameter with the crown-like morphology of the virion as observed using electron microscopy [12, 13]. The first strain of SARS coronavirus species was identified in 2002–3 [14]. As a spillover disease, it has spread throughout the world through a variety of hosts, including bats, civets, camels, and among others [15, 16]. The huge subfamily of CoVs belongs to the Nidovirales family, divided further into four primary genuses: alpha, beta, gamma and delta-coronaviruses, among which SARS CoV-2 belongs to β-coronaviruses [4, 17]. Patients with the condition may experience a variety of clinical manifestations, ranging from minor to severe symptoms such as headaches, shortness of breath, muscular pains, fever and tiredness, to multi-organ failure [18, 19]. Inhalation of respiratory droplets, usage of personal belongings of the infected patients, direct or indirect physical exposure and contamination have been regarded as the ways for the illness to spread among human beings [20].

Despite the devastating effects, no viable medicines to tackle COVID-19 have been identified so far. There are currently sixty-three diversified vaccines containing inactivated viruses (Covaxin), adeno-based vaccines (Gam-Covid by Gameleya Institute), recombinant (Novavax), mRNA (Moderna) and live attenuated vaccines that have been approved on a global basis by regulatory authorities in order to effectively diminish or uproot the effect of COVID-19. More than twenty vaccinations are in phase III clinical investigations and will be introduced with better efficacy and safety in the future. However, repurposing of existing anti-viral drugs and development of vaccines could not act straightforwardly proving themselves an effective therapy to prevent the virus for making its way to animal species through binding with spike glycoprotein. Additionally, they bring noticeable ill effects on the imperative physiological processes of human body upon administration such as revering to the deficit of potential therapy in this line of approach; there is an urgent need for the development of promising candidates [21, 22].

To become an FDA-approved drug, the molecules must have to pass through the hectic discovery process through conventional methods at the cost of billions of US dollars that too might not guarantee a potential therapy. On the other hand, drug discovery through the cost effective virtual screenings could serve the purpose of novel drug discovery for the eradication of SARS CoV-2 proving itself a need of an hour [22-25]. In order to get an effective therapeutic approach, critical insights into the viral entry and replication cycle into host’s body have been urged to understand the pathogenesis of SARS CoV-2 (Fig. 1). The viral replication process starts with the entry of virion into host cells via spike glycoprotein (S), which is followed by the use of host cell machinery to assemble and reproduce multiple viral parasites. Spike protein, which is made up of two essential subunits called S1 and S2, facilitates in virus binding and fusion to the host cell membrane, which is accomplished via endocytosis or cleavage at the interface of the two subunits by host proteases such as transmembrane protease serine type 2 (TMPRSS2) and cathepsin. Furthermore, the released 30 kb positive-sense single-stranded ssRNA (5′-3′ UTR) into the host cytoplasm uses host cell machinery to translate into polyproteins, which are then proteolytically cleaved into sixteen distinct non-structural proteins by the proteases encoded by the virus (NSPs). These NSPs construct a replicase–transcriptase complex, which incorporates RNA-dependent RNA polymerase (RdRp), helicases, endo and exonucleases, and other important enzymes involved in nucleic acid metabolism. Open reading frames required for transcription into structural proteins such as spike (S), membrane (M), nucleocapsid (N) and envelope (E) are encoded at the 3′-end of the genome (E). The proteins generated in the endoplasmic reticulum (ER) of the host cell are transformed into new viral offspring, which are then released into ER-golgi compartments and infect additional host cells [26, 27] .

Fig. 1

Life cycle of SARS CoV-2 from entry to replication into host cells

In search of SARS CoV-2 inhibitors, Buchwald and co-workers have performed the virtual screening of library consisting of organic dyes and reported p-nitrobenzamide derivative (A, Fig. 2) with the IC50 of 5.6 µM, as spike-ACE2 protein–protein interaction blockers inhibiting the attachment and entry of coronavirus to human cells [28]. The multidisciplinary approach performed by Haselhorst et al. via screening of 57,641 ligands to inhibit human ACE2 or spike protein led to identification of evans blue (B, Fig. 2) as potent anti-viral agent with IC50 of 28.1 µM against vero-E6 cells incubated with SARS CoV-2 [29]. Three different peptide fragments composed of N-terminal amino acid residues of α1 helix of ACE2 peptidase domain have been reported by de Olivera and co-workers as effective agents to inhibit the binding of SARS CoV-2 spike protein with human ACE2 as identified through molecular modelling [30]. Further, Baig et al. identified the 13-amino acid peptide (FLDKFNHEAEDLF) as spike protein inhibitor with minimal inhibition (40% at 100 µM concentration) of ACE2-spike protein–protein interactions through enzyme-linked immunosorbent (ELISA)-dependent inhibitory assay [31].

Fig. 2

Reported spike-ACE2 inhibitors inhibiting SARS CoV-2 attachments and their entries in host cells

The use of computational studies has promoted the design and discovery of new chemicals scaffolds acting against the choice of drug targets through molecular modelling techniques like molecular docking, ADMET analysis and molecular dynamics (MD) simulation [32]. Recently, we have reported the virtual screening of small molecules against RdRp [33] and Mpro [34] using molecular docking, ADMET analysis and MD simulations. In continuation towards our endeavour for search of SARS CoV-2 inhibitors, we become interested for the in silico-based virtual screening of 175,851 ligands of Asinex BioDesign library considering the massive demand for SARS CoV-2 inhibitors having good binding affinity towards the spike glycoprotein.

Results and discussion

Pharmacophore hypothesis

SiteMap analysis

The protein structure of spike protein (PDB ID: 6VXX) having resolution of 2.8 Å with the symmetrical trimeric structure possessing two receptor binding domains (S1 and S2 domains) has been deposited in protein data bank by Veesler et al. [35]. There has not been any ligand co-crystallized with the protein of interest. As a result, attempts have been made to determine the probable binding site of spike glycoprotein using the SiteMap tool of Schrödinger [36]. Five possible active sites have been predicted, and their drugability scores (Dscore) and SiteScores have been calculated by consideration of the pocket volume, enclosure, and the degree of hydrophobicity (Table 1) [37]. Next, the site having the highest Dscore (1.103642) has been utilized for the pharmacophore hypothesis development which includes residues Lys41, Ile197, Asp198, Tyr200, Lys202, Asp228, Arg355, Cys379, Tyr380, Gly381, Val382, Ser383, Tyr396, Ala411, Pro412, Gly413, Gln414, Thr415, Ala419, Asp420, Lys424, Pro426, Asp427, Asp428, Phe429, Thr430, Phe464, Ser514, Phe515, Glu516, Leu517, Leu518, His519, Glu748, Asn751, Leu752, Gln755, Tyr756, Gln762, Leu763, Arg765, Ala766, Gly769, Ile770, Glu773, Val951, Gln954, Gln957, Ala958, Thr961, Leu962, Gln965, Phe970, Gly971, Ile973, Ser974, Asp979, Leu981, Arg983, Leu984, Asp985, Pro986, Pro987, Glu988, Ala989, Glu990, Val991, Gln992, Asp994, Arg995, Thr998, Gly999, Leu1001, Gln1002, Ser1003, Gln1005, Thr1006, Tyr1007, Val1008, Thr1009, Gln1010, Gln1011, Leu1012, Ile1013, Arg1014, Ala1015, Ala1016, Glu1017, Ile1018 and Arg1019 (Table 1, Entry 1).

Table 1

Predicted five sites from the sitemap analysis for the selected spike protein (PDB ID: 6VXX)

Entry	SiteScore	Dscore
1	1.069965	1.103642
2	1.054111	1.057707
3	1.035711	1.025884
4	1.128571	1.011641
5	1.027926	0.997882

Pharmacophore development and ligand screening

Next, using the phase module of Schrödinger, an energy-based pharmacophore (e-pharmacophore) hypothesis model has been generated around the predicted binding site of spike glycoprotein to obtain the steric and electrostatic attributes [38]. The model with seven features (RRRRRRH) has been generated comprising of six aromatic rings (R) and a hydrophobic feature (H, see supporting information, Fig. S2). Total 175,851 numbers of compounds have been obtained from Asinex BioDesign library 2020.1 from online available sources [39]. With a view of validation of the developed hypothesis, all the compounds have been processed through preparation of ligands using LigPrep [40] followed by screening using phase [38]. Total 38,267 molecules from the selected database have matched with the minimum four features of the developed hypothesis, and they have been further utilized for the molecular docking.

Molecular docking

Molecular modelling using docking [41-44] has been considered as an important computational tool to predict the binding interactions of the ligands with the active site of protein in search of their mode of action against SARS CoV-2 [45-52]. In this context, we have performed the molecular docking in the search of potent SARS CoV-2 inhibitors using the obtained hypothesis model. The receptor grid box with 10 Å has been generated to perform the molecular docking on the selected active site from the SiteMap analysis using the Grid-based Ligand Docking with Energetics (GLIDE) module of Schrödinger [53]. To get the detailed knowledge about the binding strength and types of interactions of the ligand with receptor, we have performed sequential docking at three different precision levels such as high throughput virtual screening (HTVS), standard precision (SP) and extra precision (XP). First of all, the HTVS has been performed using 38,267 ligands, qualified in phase screening, from which 1094 molecules have been found with the docking score of ≤ − 4.6. Next, they have been further processed for docking with SP mode and 175 molecules have been identified with the docking score ≤ − 5.9 in SP docking which have been passed through the XP docking with the generation of 10 poses for each ligand molecule. From that, seventeen compounds with highest docking scores (≤ − 6.0) have been identified and summarized in Table 2, and 2-D interactions are presented in Figs. 3, 4 and 5.

Table 2

Molecular docking result of identified hits with the docking score ≤ 6 in XP

Fig. 3

2D interactions of compounds having docking score of ≤ − 6.0 using XP module against spike glycoprotein. HB has been represented as purple-coloured arrows and π–π bond as green lines; solvent exposures have been presented through grey spot and salt bridge as red-blue lines

Fig. 4

2D interactions of compounds having docking score of ≤ − 6.0 using XP module against spike glycoprotein. HB has been represented as purple-coloured arrows and π–π bond as green lines; solvent exposures have been presented through grey spot and salt bridge as red-blue lines

Fig. 5

2D interactions of compounds having docking score of ≤ − 6.0 using XP module against spike glycoprotein. HB has been represented as purple-coloured arrows and π-π bond as green lines; solvent exposures have been presented through grey spot and salt bridge as red-blue lines

Further, we have studied the detailed 3D interactions of the docked ligands having the docking score ≤ − 6.0 using the PyMol 2.4.0 [54]. Compound 1 has shown the highest docking score of − 7.34 by formation of the hydrogen bond (HB) interactions with Phe970 and Asp994 along with the two salt bridges with Asp994 as presented in Fig. 6a. Hydroxy group and nitrogen atoms were involved in the bond formation. The nitrogen atom present in the five membered ring of compound 2 has interacted with Arg995 through π-cation interaction and ionic interactions (salt bridge) with Asp994 (Fig. 6b). Another compound 3 has also shown the similar interactions with Asp994 residue via formation of salt bridge with nitrogen atom (Fig. 6c) being tight binding in the predicted active site of spike glycoprotein. Benzene rings of compound 4 (Fig. 6d) and nitrogen atom present in the benzene ring of 5 (Fig. 6e) have been found to form double and single and π-π interactions, respectively, with Tyr756. Additionally, compound 6 has been fitted in the active site with maximum exposures to solvent molecules specifically with the fluorine atom (Fig. 6f).

Fig. 6

3D interactions of identified hits including 1 (a), 2 (b), 3 (c), 4 (d), 5 (e) and 6 (f) through XP docking. The poses for the represented ligands (yellow coloured ball and stick models) and protein (coloured cartoons) have been generated and visualized using PyMol 2.4.1 [54]

Hydroxy and amine group of 1H-imidazole derivative (7) has been involved in π-cation interaction with Arg995 along with the formation of salt bridge (Asp994) and HB (Asp994 and Phe970, Fig. 7a). Biologically divergent 1,2,4-oxadiazole derivative (8) [55] has been tightly accommodated with several solvent molecules and surrounded by amino acid residues in the active site of spike glycoprotein (Fig. 7b). The nitrogen atoms of hybrid heterocyclic scaffold 9 have been found to form HB with Thr998 (Fig. 7c) with docking score of − 6.53 (Entry 9, Table 2). Compound 10 has been found with formation of HB with Thr998 and Thr970, salt bridge with Asp994, π-cat with Arg995 and π–π bond with Tyr756 (Fig. 7d) showing the string interactions in the active site of protein. Hydroxy and amine groups were found to involve in the interaction. Pyrimidine derivative compound 11 has been found to fit into the active site by formation of HB with Thr998 along with two π-π bonds with Tyr756 (Fig. 7e). The HB formation has been observed between Gln1002 and ethereal oxygen (HB acceptor) attached to 4-chlorophenyl ring of compound 12 (Fig. 7f) having docking score of -6.33 (Entry 12, Table 2).

Fig. 7

3D interactions of identified hits including 7 (a), 8 (b), 9 (c), 10 (d), 11 (e) and 12 (f) through XP docking. The poses for the represented ligands (yellow coloured ball and stick models) and protein (coloured cartoons) have been generated and visualized using PyMol 2.4.1 [54]

Compound 13 has been found to form the π-cation interaction with Arg995 (Fig. 8a) along with formation of salt bridge and HB formation with Asp994 resulting in docking score of − 6.23 (Entry 13, Table 2). Nitrogen of the amine group of non-heterocyclic compound 14 has formed HB with Thr998 residue and π–π bond interaction with Tyr756 residue (Fig. 8b), whereas nitrogen atom and hydroxy group present in compound 15 have been found to forms HB with Phe970 and Arg995 residues (Fig. 8c). Synthetically divergent indole [56] compound 16 has been found in the active site with mainly solvent exposure (Fig. 8d) with the docking score of − 6.04 (Entry 16, Table 2). Nitrogen of the amine group present in the 2-ethyl imidazole derivative 17 has been observed with the formation of HB with Thr998 and π-π bonds with the Tyr756 residues of both chains A and B (Fig. 8e).

Fig. 8

3D interactions of identified hits including 13 (a), 14 (b), 15 (c), 16 (d) and 17 (e), through XP docking. The poses for the represented ligands (yellow coloured ball and stick models) and protein (coloured cartoons) have been generated and visualized using PyMol 2.4.1 [54]

ADMET analysis

Further, we have performed the analysis of absorption, distribution, metabolism, excretion and toxicity (ADMET) parameters for the identified hits through XP docking using QikProp module of Schrödinger [57]. Various physicochemical properties like molecular weight (MW), oil/water partition coefficient (LogP), water solubility (LogS), IC50 value for blockage of (human ether-a-go-go-related gene (HERG) K+ channels, gastrointestinal (GI) cell for cell permeability (QPP caco2), brain/blood partition coefficient (Log BB), number of metabolic reactions (Metab), binding affinity to human serum albumin (QPLog Khsa) and human oral absorption in percentage (% HOA) have been studied for the selected compounds and are presented in Table 3. Mostly all the hits have satisfied the ADMET criteria, and all the hit molecules may have good oral absorption. This provides better scope for the drug molecule development in the future.

Table 3

ADMET parameters of the identified hits

Compd. No	MWa	LogPb	Log Sc	QPlog HERGd	QPP caco2e	Log BBf	Metabg	QPLog Khsah	% HOAi
1	453.56	5.18	− 3.90	− 5.69	591.32	− 0.16	6	0.83	93.95
2	483.58	5.96	− 6.25	− 7.59	760.69	− 0.26	7	1.03	100
3	362.49	4.35	− 5.49	− 7.58	401.76	− 0.13	4	0.76	100
4	449.03	5.06	− 3.73	− 8.30	287.78	0.54	6	0.94	87.64
5	397.44	5.34	− 5.57	− 7.05	1267.41	0.70	6	0.85	100
6	391.37	4.57	− 5.60	− 7.03	746.09	0.66	4	0.60	100
7	453.56	4.77	− 3.58	− 6.41	685.05	− 0.18	5	0.61	100
8	373.38	3.88	− 4.46	− 6.58	524.52	0.43	4	0.41	100
9	414.51	3.61	− 4.30	− 5.93	194.37	− 0.31	3	0.47	89.06
10	451.56	5.43	− 4.82	− 6.42	682.53	− 0.28	7	0.97	96.48
11	407.51	4.59	− 4.54	− 7.04	716.06	− 0.09	6	0.68	100
12	484.04	6.13	− 6.27	− 7.25	918.95	− 0.07	7	1.14	100
13	470.01	5.63	− 4.56	− 6.12	1244.17	0.23	5	0.80	100
14	407.48	5.41	− 5.51	− 7.90	959.10	0.14	3	0.91	100
15	469.56	5.21	− 4.48	− 6.20	713.70	− 0.16	7	0.82	95.57
16	416.49	6.73	− 8.63	− 7.49	3518.97	− 0.27	6	1.41	100
17	467.58	5.18	− 5.38	− 6.81	872.36	− 0.13	6	1.06	100
Standard value	< 500	2.0–6.5	> − 4	> − 5	< 25 (poor) and > 500 (excellent)	–	1–8	– 1.5 to 1.5	0–100%

aMW: Molecular weight; (Da); bOil/water partition coefficient; cSolubility; dIC50 value for blockage of HERG K+ channels; eGut–blood barrier in nm/s permeability of the cell model; fBrain/blood partition coefficient; gNumber of metabolic reactions; hBinding to human serum albumin; iHuman oral absorption. The reference values for ADMET have been taken from the Ref [58]

The Lipinski’s rule of five (Ro5) has been recognized as the suitable in silico tool for the assessment of essential physico-chemical properties of the hits or leads identified during pre-clinical settings [59]. To ensure the drug likeliness of the hit molecules, the physico-chemical properties for Ro5 have been also evaluated using the QikProp module of Schrödinger [57]. The parameters described under Lipinski’s rule of five such as molecular weight (MW), number of hydrogen bond donor (HBD) and acceptor (HBA), oil to water partition coefficient (LogP) and number of rotatable bonds (RB) have been studied (Table 4). The analysis of these properties revealed that none of the hits has violated more than one rule as per the criteria of the Lipinski’s Ro5, highlighting the significance of these hits in the search of SARS CoV-2 inhibitors.

Table 4

Analysis of physico-chemical properties using Lipinski’s rule of Five (Ro5)

Compd. No	MWa	HBAb	HBDc	LogPd	RBe	Violationf
1	453.56	5.75	1	5.18	11	1
2	483.58	6.5	1	5.96	12	1
3	362.49	4.5	2	4.35	5	0
4	449.03	5.2	2	5.06	9	1
5	397.44	4.5	0	5.34	4	1
6	391.37	5	0	4.57	2	0
7	453.56	6.25	1	4.77	11	0
8	373.38	5	0	3.88	2	0
9	414.51	6	2	3.61	5	0
10	451.56	6.5	1	5.43	11	1
11	407.51	6.5	1	4.59	8	0
12	484.04	5.75	1	6.13	11	1
13	470.01	5.75	0	5.63	10	1
14	407.48	4.45	2	5.41	8	1
15	469.56	6.5	1	5.21	11	1
16	416.49	3.25	2	6.73	6	1
17	484.04	5.75	1	6.13	11	1
Standard value	< 500	< 10	< 5	< 5	< 5	–

aMolecular weight (Da), bNumber of hydrogen bond acceptors, cNumber of hydrogen bond donors, dPartition coefficient in oil to water, eNumber of rotational bonds, fNumber of violated parameters from Ro5

The Brain Or IntestinaL EstimateD permeation method (BOILED-Egg) gives the estimation of accessibility of compounds to gastrointestinal (GI) tract and blood–brain barrier [60]. So, we also evaluated the accessibility of the hits using SwissADME [61]. The boiled egg model revealed all the hit molecules possessed satisfactory GI absorption along with inhibition of the P-glycoprotein, a protein responsible for efflux of drugs from cells (Fig. 9). All the compounds (except 3 and 16) may have sufficient permeability across the blood–brain permeability (BBB) indicating their usefulness to treat the infectious disorders related to central nervous system.

Fig. 9

The BOILED-egg model of identified hit molecules (1–17) obtained through SwissADME. BBB and GI permeability have been indicated through yellow and colourless regions, respectively. The blue circles denote inhibition of P-glycoprotein by the corresponding hits

The presently identified compounds (1–17) gain the advantage of sufficient metabolic stability and free from nitro-aryl derived toxicity over the literature reported compounds (A and B, Fig. 2), whereas these organic (azo)dyes (A and B) suffer from the limitations of intense colour and metabolic instabilities limiting their therapeutic potential against pathogenic diseases [62]. Further, compound A possesses nitro-aryl features imparting the potential risk of mutagenicity due to interaction of nitro group with DNA [63]. Peptides identified by de Oliveira et al. [30] and Baig et al. [31] may not have sufficient oral bioavailability as compared to identified small molecules (1–17) due to degradation of these peptides under the acidic environment of stomach and peptidase enzymes [64].

MD simulation

Molecular dynamics simulations have been recognized as significant tool to claim the stability of the ligand and to reveal macromolecular structure to function relationships in the field of drug discovery [65-68]. Henceforth, we performed the molecular dynamics (MD) simulations at various time points up to 10 ns using GROMACS 2020.1 to assess the stability of the hits into the complexes obtained from ADMET and Lipinski rules analysis [69, 70]. The graphical representation of plots of statistical parameters obtained by the MD simulation of compound 1 is presented in Fig. 10. The ligand–receptor complex of 1 with spike was found with root mean square deviation (RMSD) value ranging from 5.81 to 10.14 nm with an average of 8.889 nm (Fig. 10a) for the ligand and 7.06–10.14 nm with an average of 9.710 nm (Fig. 10b) for the protein. Although these values were little high, the stability of the ligand into the active site was found consistent as evident from Fig. 10a, b. The radius of gyration (RoG) for the same complex ranging from 7.13 to 8.90 nm with an average of 8.311 nm (Fig. 10c) indicated the compactness of the complex. The surface area accessed by the solvent molecules was found within the range from 1550 to 1590 nm2 with an average of 1586.597 nm2 (Fig. 10d). The plot of number of hydrogen bonds (HB) vs simulation time revealed maximum five HBs between the ligand and receptor within the time period of 10 ns (Fig. 10e). The electrostatic (coulombic short-range, Coul-SR) and van der Waals/hydrophobic interactions (LJ-SR) for the complex were found − 230.361 ± 7.5 and − 146.553 ± 5 kJ/mol, respectively, which indicated the key role of the electrostatic interactions to stabilize the complex of 1 with spike glycoprotein.

Fig. 10

The schematic plots of (6a) RMSD-L, (6b) RMSD-P, (6c) RoG, (6d) SASA and (6e) HB for the complex of compound 1 with spike glycoprotein

The graphical representation of plots of statistical parameters obtained by the MD simulation of compound 2 is represented in Fig. 11. The ligand–receptor complex of 2 with the spike glycoprotein was found with RMSD value ranging from 7.08 to 11.32 nm with an average of 9.05 nm (Fig. 11a) for the ligand and 7.39–7.80 nm with an average of 7.70 nm (Fig. 11b) for the protein. Although these values were high, the stability of the ligand into the active site was found consistent during 0–8 ns, and little instability was observed after 8 ns in the RMSD. The radius of gyration (RoG) for the same complex ranging from 8.4 to 8.8 nm with an average of 8.586 nm (Fig. 11c) indicated the compactness of the complex. The surface area accessed by the solvent molecules was found within the range from 1030 to 1070 nm2 with an average of 1060.85 nm2 (Fig. 11d). The plot of number of hydrogen bonds (HB) vs simulation time revealed maximum four HBs between the ligand and receptor within the time period of 10 ns (Fig. 11e). The electrostatic (coulombic short-range, Coul-SR) and van der Waals/hydrophobic interactions (LJ-SR) for the complex were found − 87.9052 ± 18 and − 72.2626 ± 13 kJ/mol, respectively.

Fig. 11

The schematic plots of (7a) RMSD-L, (7b) RMSD-P, (7c) RoG, (7d) SASA and (7e) HB for the complex of compound 2 with spike glycoprotein

Conclusions

In the search of effective SARS CoV-2 inhibitors inhibiting the viral entry of coronavirus into host cells, the systematic virtual screening of ligands from Asinex BioDesign library 2020.1 has been performed against spike glycoprotein, a promising target of SARS CoV-2. The predicted binding site of spike glycoprotein has been used for pharmacophore-based screening, wherein the selected ligands with essential pharmacophoric features have been subjected for three different levels of precise molecular docking to shortlist the identified hits with the desired docking score. Total seventeen hits having the promising docking score have qualified the ADMET parameters and Lipinski’s rule of five properties to claim their oral bioavailability. The dynamics simulation studies have demonstrated the stability of identified hits (1 and 2) into the predicted catalytic site of protein with sufficient compactness, uniqueness and stability. In summary, the present studies accomplishes that the identified hits may interfere with the key residues of spike glycoprotein effectively as SARS CoV-2 inhibitors with sufficient binding affinity, good stability into the active site with the acceptable drug likeliness. These studies have provided SARS CoV-2 inhibitors having the potential of inhibiting the viral entry into the host cell by acting on the viral spike protein against COVID-19.

Computational details

Preparation of ligands

Asinex BioDesign library 2020.1, comprising of 175,851 molecules with pharmacologically important structural features in the chemical scaffolds with synthetic feasibility, has been accessed through open source databases [39]. The ligands were prepared for molecular modelling using LigPrep, where the possible states of ionization were generated using Epik at physiological pH (7.0 ± 2.0) with consideration of metal binding states, removal of salt forms and generation of tautomeric isomers or stereoisomers using the OPLS3e force field [40].

Preparation of proteins

The 3D structure of SARS CoV-2 spike glycoprotein (PDB ID: 6VXX) having resolution of 2.8 Å has been accessed from RCSB protein data bank [71]. The protein has been modified using protein preparation wizard (PrepWizard) of Schrödinger suite (Maestro 12.5) [72]. The hydrogen consistency, steric relations, bond orders and total charges have been generated followed by optimization and energy minimization of protein using the force field (OPLS3e) to assure the structural accuracy of final protein.

Site map analysis

The selected protein does not contain any co-crystallized ligand or inhibitor. In a view of this, the site map analysis has been performed to identify specific receptor binding site [36]. The possible binding site regions with minimum fifteen site points have been generated under more restrictive hydrophobic environment using standard grid to visualize and evaluate the top-most binding sites. For further molecular modelling, binding site with the highest drugability score (DScore) has been selected as an active site from the nominated top five binding sites.

Pharmacophore development and screening

The selected protein of interest (PDB ID: 6VXX) does not contain any co-crystallized ligand [35]. The residues have been specified through the atom specification language (ASL) for the selected best active site with highest D score followed by the generation of receptor cavity-based E-pharmacophore hypothesis using the phase module of Schrödinger [38]. The model has been generated for the selected active site along with the generation of seven features (RRRRRRH), comprising of aromatic ring (R) and hydrophobic feature (H). The molecules prepared after the LigPrep have been subjected to phase ligand screening to fit with minimum four of the identified features of the generated pharmacophore hypothesis model. Maximum fifty conformers for each ligand have been generated and subjected to energy minimization. The phase screen score for each conformer has been calculated based on preset acceptor and donor as negative and positive equivalent, respectively. Total 38,267 molecules have qualified with the set features of hypothesis from the selected drug library.

Receptor grid generation

The grid box of size 10 Å has been generated using GLIDE module of Schrödinger on the identified binding site from site map analysis with the highest DScore [53]. The van der Waals radius scaling factor and the partial charge cut-off have been set 1 and 0.25, respectively, to soften the potential for nonpolar parts of the protein.

Ligand docking

Next, molecular docking has been performed using GLIDE module of Schrödinger. The 38,267 ligands which have qualified through the phase screening, have been docked on selected active sites of the protein using high throughput virtual screening (HTVS) through flexible ligand sampling and one pose per each ligand has been generated. A total of 1094 molecules with docking score of ≤ − 4.6 obtained in HTVS docking have been further docked on the same active site using standard precision (SP) to generate one pose for each ligand. Again 175 top compounds having docking score ≤ − 5.9 have been further docked with extra precision (XP) to generate 10 poses per ligand. ADMET analysis with Lipinski parameters has been performed for 17 hits having docking score of ≤ − 6.0 in XP docking. The interactions of the ligand with the amino acids of the active site have been visualized using PyMol 2.4.1 [54].

ADMET analysis and Lipinski rule of five

The absorption, distribution, metabolism, excretion and toxicity (ADMET) analysis and assessment of Lipinski parameters have been carried out using QikProp module of Schrödinger for the obtained 16 hits which have achieved the docking score less than − 6.33 in XP docking [57]. The various physical, chemical and functional group properties such as molecular weight, oil/water partition coefficient, solubility, IC50 value for blockage of HERG K+ channels, gut–blood barrier in nm/s permeability of the cell model, brain/blood partition coefficient, number of metabolic reactions, binding to human serum albumin and human oral absorption have been studied from the manual analysis of the result of ADMET analysis. The drug likeliness features like molecular weight, hydrogen bond donors, hydrogen bond acceptors, partition coefficient (oil/water) and number of rotatable bonds have been studied through Lipinski rule of five.

Molecular dynamics (MD) simulation

Hit molecule obtained from the manual analysis of ADMET and Lipinski parameters have been incorporated to perform the MD simulation using GROningen MAchine for Chemical Simulations (GROMACS) 2020.1 [69, 70] software. CHARMM36 (Chemistry at Harvard Macromolecular Mechanics) was used as an all atom force field [73] and CHARMM General Force Field (CGenFF) server [74, 75] was used to generate the topology of the ligands. After solvation (TIP3P water model) and neutralization (sodium chloride), the complex was subjected to energy minimization followed by equilibration through the canonical isovolumetric-isothermal (NVT) and isobaric-isothermic (NPT) canonicals for 100 ps. The energy minimized complexes have been subjected for final MD run of 10 ns for successful completion.

Below is the link to the electronic supplementary material.

Supplementary file1 (PDF 419 KB)