key: cord-0882097-k7a4ui3i
authors: Dhameliya, Tejas M.; Nagar, Prinsa R.; Gajjar, Normi D.
title: Systematic virtual screening in search of SARS CoV-2 inhibitors against spike glycoprotein: pharmacophore screening, molecular docking, ADMET analysis and MD simulations
date: 2022-02-08
journal: Mol Divers
DOI: 10.1007/s11030-022-10394-9
sha: 5a7f036538fa90f8b085d0887e0315e4da48f5ef
doc_id: 882097
cord_uid: k7a4ui3i

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. GRAPHICAL ABSTRACT: [Image: see text] SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s11030-022-10394-9.

Human ether-a-go-go-related gene HOA Human oral absorption HTVS High-throughput virtual screening M Membrane protein N Nucleocapsid nCoV- 19 Novel coronaviruses NSPs Non-structural proteins OPLS3e

Optimized potentials for liquid simulations pp1a/b Polyproteins RdRp RNA-dependent RNA polymerase S Spike glycoprotein

The novel coronavirus disease 2019 , 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] [4] [5] [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] [23] [24] [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 ERgolgi compartments and infect additional host cells [26, 27] .

In search of SARS CoV-2 inhibitors, Buchwald and coworkers have performed the virtual screening of library consisting of organic dyes and reported p-nitrobenzamide derivative (A, Fig. 2 ) with the IC 50 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 IC 50 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 ACE2spike protein-protein interactions through enzyme-linked immunosorbent (ELISA)-dependent inhibitory assay [31] .

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 

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).

Next, using the phase module of Schrödinger, an energybased 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 (RRR RRR H) 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 modelling using docking [41] [42] [43] [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] [46] [47] [48] [49] [50] [51] [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.

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 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 π-π 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) .

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 ). 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) .

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), IC 50 value [54] 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.

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. 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.

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 

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] [66] [67] [68] . Henceforth, Fig. 9 The BOILED-egg model of identified hit molecules (1-17) obtained through Swis-sADME. 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 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 nm 2 with an average of 1586.597 nm 2 (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 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 electrostatic interactions to stabilize the complex of 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 nm 2 with an average of 1060.85 nm 2 (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.

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 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 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 pharmacophorebased 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.

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 Lig-Prep, 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] .

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.

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.

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 (RRR RRR H), 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.

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.

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] .

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, IC 50 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.

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.

The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s11030-022-10394-9.

Characteristics of SARS-CoV-2 and COVID-19

Origin and evolution of pathogenic coronaviruses

Insights into the recent 2019 novel coronavirus (SARS-CoV-2) in light of past human coronavirus outbreaks

The recent outbreaks of human coronaviruses: a medicinal chemistry perspective

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modeling study

The continuing 2019-nCoV epidemic threat of novel coronaviruses to global health: the latest 2019 novel coronavirus outbreak in Wuhan, China

Coronavirus disease 2019 (COVID-19)

WHO Coronavirus Disease

The socio-economic implications of the coronavirus pandemic (COVID-19): a review

A brief review of socio-economic and environmental impact of COVID-19

Genotype and phenotype of COVID-19: their roles in pathogenesis

Current perspective of antiviral strategies against COVID-19

Coronaviruses. In: Baron S (ed) Medical Microbiology. University of Texas Medical Branch at Galveston

A novel coronavirus from patients with pneumonia in China

Characterization and noncovalent inhibition of the deubiquitinase and deISGylase activity of SARS-CoV-2 papain-like protease

Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding

Coronavirus diversity, phylogeny and interspecies jumping

Clinical characteristics of coronavirus disease 2019 in China

The epidemiology and pathogenesis of coronavirus disease (COVID-19) outbreak

A review of coronavirus disease-2019

Drug development and medicinal chemistry efforts toward SARS-coronavirus and COVID-19 therapeutics

Computational drug discovery and repurposing for the treatment of COVID-19: a systematic review

COVID-19 drug repurposing: A review of computational screening methods, clinical trials, and protein interaction assays

A review on potential of natural products in the management of COVID-19

The race to treat COVID-19: Potential therapeutic agents for the prevention and treatment of SARS-CoV-2

Coronavirus membrane fusion mechanism offers a potential target for antiviral development

Host-pathogen interaction in COVID-19: pathogenesis, potential therapeutics and vaccination strategies

Small-molecule inhibitors of the coronavirus spike: ACE2 protein-protein interaction as blockers of viral attachment and entry for SARS-CoV-2

Multidisciplinary approaches identify compounds that bind to human ACE2 or SARS-CoV-2 spike protein as candidates to block SARS-CoV-2-ACE2 receptor interactions

Shedding light on the inhibitory mechanisms of SARS-CoV-1/CoV-2 spike proteins by ACE2-designed peptides

A novel therapeutic peptide blocks SARS-CoV-2 spike protein binding with host cell ACE2 receptor

Computer-aided drug design. In: Poduri R (ed) Drug discovery and development

In search of SARS CoV-2 replication inhibitors: virtual screening, molecular dynamics simulations and ADMET analysis

In search of RdRp and Mpro inhibitors against SARS CoV-2: molecular docking, molecular dynamic simulations and ADMET analysis

Structure, function, and antigenicity of the SARS-CoV-2 spike glycoprotein

Schrödinger Release 2020-3: SiteMap, Schrödinger, LLC

Druggability analysis and structural classification of bromodomain acetyllysine binding sites

Schrödinger Release 2020-3: Phase, Schrödinger, LLC

Schrödinger Release 2020-3: LigPrep, Schrödinger, LLC

Benzo[d]thiazole-2-carbanilides as new anti-TB chemotypes: design, synthesis, biological evaluation, and structure-activity relationship

Synthesis, biological evaluation and structure-activity relationship of 2-styrylquinazolones as anti-tubercular agents

N-Arylalkylbenzo[d] thiazole-2-carboxamides as anti-mycobacterial agents: design, new methods of synthesis and biological evaluation

Identification of anti-mycobacterial agents against mmpL3: virtual screening, ADMET analysis and MD simulations

Quinolines-based SARS-CoV-2 3CLpro and RdRp inhibitors and Spike-RBD-ACE2 inhibitor for drug-repurposing against COVID-19: an in silico analysis

Ligand and structure-based virtual screening applied to the SARS-CoV-2 main protease: an in silico repurposing study

Ultrasound assisted synthesis of 3-alkynyl substituted 2-chloroquinoxaline derivatives: their in silico assessment as potential ligands for N-protein of SARS-CoV-2

High throughput virtual screening to discover inhibitors of the main protease of the coronavirus SARS-CoV-2

Computational determination of potential inhibitors of SARS-CoV-2 main protease

Investigation of some antiviral N-heterocycles as COVID 19 drug: molecular docking and DFT calculations

Putative inhibitors of SARS-COV-2 main protease from a library of marine natural products: a virtual screening and molecular modeling study

Pyrrolo[2,3-b] quinoxalines in attenuating cytokine storm in COVID-19: their sonochemical synthesis and in silico/in vitro assessment

The PyMOL molecular genetics graphics system, DeLano Scientific LLC, San Carlos 55

Synthetic account of indoles in search of potential anti-mycobacterial agents: a review and future insights

Schrödinger Release 2020-3: QikProp, Schrödinger, LLC

Structure-based design, synthesis and biological evaluation of a newer series of pyrazolo[1,5-a]pyrimidine analogues as potential anti-tubercular agents

Experimental and computational approaches to estimate solubility and permeability in drug discovery and devlopment settings

A BOILED-Egg to predict gastrointestinal absorption and brain penetration of small molecules

SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules

Metabolism of aZO dyes: implication for detoxication and activation

Mutagenicity of nitroaromatic compounds

Advancements in docking and molecular dynamics simulations towards ligand-receptor interactions and structure-function relationships

Molecular dynamics simulations in drug discovery and pharmaceutical development

Use of molecular dynamics simulations in structure-based drug discovery

Accelerating COVID-19 research using molecular dynamics simulation

GROMACS 2020.1 (Manual Version 2020.1) Zenodo

GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers

Protein Preparation Wizard

CHARMM36m: an improved force field for folded and intrinsically disordered proteins

CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields

Extension of the CHARMM general force field to sulfonyl-containing compounds and its utility in biomolecular simulations

The authors would like to acknowledge the support provided by Schrödinger Inc. for performing molecular modelling under the banner of Drug Discovery Hackathon 2020 (DDH-2020), a joint initiative of All India Council of Technical Education (AICTE), New Delhi, Council of Scientific and Industrial Research (CSIR), New Delhi, and Government of India.Author contributions All the authors contributed equally for this work. TMD performed docking, data generation, conceptualization, supervision, writing-original draft and revised manuscript preparation, and reviewing and editing; PRN contributed to methodology, data generation, docking, writing-original draft and revised manuscript preparation; and NDG contributed to methodology, data generation, docking, ADMET assay, MD simulations, writing-original draft and revised manuscript preparation.