Imidazoles and benzimidazoles as putative inhibitors of SARS-CoV-2 B.1.1.7 (Alpha) and P.1 (Gamma) variant spike glycoproteins: A computational approach.

COVID-19 is an unprecedented pandemic threatening global health, and variants were discovered rapidly after the pandemic. The two variants, namely the SARS-CoV-2 B.1.1.7 (Alpha) and P.1 (Gamma), were formed by the mutations in the receptor binding domain of spike glycoprotein (SGP). These two variants are known to possess a high binding affinity with the angiotensin-converting enzyme 2. Amidst the rapid spread of these mutant strains, research and development of novel molecules become tedious and labour-intensive. Imidazole and benzimidazole scaffolds were selected in this study based on their unique structural features and electron-rich environment, resulting in increased affinity against a variety of therapeutic targets. In the current study, imidazole- and benzimidazole-based anti-parasitic drugs are repurposed against SARS-CoV-2 Alpha and Gamma variant spike glycoproteins using computational strategies. Out of the screened 15 molecules, flubendazole and mebendazole have exhibited promising binding features to the two receptors (PDB ID: 7NEH and 7NXC), as evidenced by their glide score and binding free energy. The results are compared with that of the two standard drugs, remdesivir and hydroxychloroquine. Flubendazole and mebendazole have become convenient treatment options against mutant lineages of SARS-CoV-2. The edge of the flubendazole was further established by its stability in MD simulation conducted for 100 ns employing GROMACS software. Further, in vitro and in vivo studies are essential to understand, if flubendazole and mebendazole indeed hold the promise to manage SARS-CoV-2 mutant stains. SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s11696-021-01900-8. © Institute of Chemistry, Slovak Academy of Sciences 2021.

Introduction

The rampant dissemination of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) or COVID-19 amongst humans hastened its molecular evolution (Hui et al. 2020). SARS-CoV-2 has been accumulating mutations at a rate of up to two nucleotide modifications per month. Current SARS-CoV-2 isolates have at least 20 nucleotide variations in their genomic sequences than the earliest isolates (Faria et al. 2021). The majority of the heterogeneity observed in SARS-CoV-2 isolates is due to mutations in the receptor binding domain (RBD) of spike glycoprotein (SGP), which facilitates virion adhesion to the angiotensin-converting enzyme 2 (ACE2) on human cell surfaces (Shang et al. 2020, Burki 2021).

The appearance of the SARS-CoV-2, B.1.1.7 lineage (UK strain) (Alpha), discovered in early December 2020 by the UK Consortium of COVID-19 Genomics, is one example of rapid viral (Kupferschmidt 2020) molecular evolution . The Alpha variant is interesting because it has 17 mutations in its genomic sequence, and 8 out of them are found in the gene that encodes the SGP. The Alpha lineage phenotype has also received recognition because of its rapid transmission rates in humans than the other recognized SARS-CoV-2 lineages. N501Y has been recognized as a residue that enhances the binding affinity with human ACE2. As per the studies by Faria et al. 2021, it is clear that the N501Y mutation observed in the Alpha lineage at the Spike gene covers one of the six primary interacting amino acid residues within the SGP (Starr et al. 2020; Sternberg and Naujokat 2020).

One more strain of SARS-CoV-2, P.1 lineage (Brazil strain) (Gamma), or 20 J/401Y.V3 variant, was first identified by the National Institute of Infectious Diseases (NIID) on 6 January 2021 in Japan (Wang et al. 2021; Maggi et al. 1249). This variant also has 17 mutations, 10 of which are in its SGP. However, three mutations were designated to be of particular concern: K417T, E484K, and N501Y. In March 2021, Gamma lineage has been designated as a deadly variant (10–80% lethal), producing nearly 10 times more viral load than in other COVID-19 patients and exhibited 1.4–2.2 times higher transmissibility. Besides, it can evade 25–61% of immunity from previous SARS-CoV-2 diseases, increasing the probability of recurring infection. These two mutant lineages were classified as variants of concern by Centers for Drug Control and Prevention (CDC). In addition, WHO has currently discovered some standard therapies using some drug combinations, investigational drugs, and vaccines (Villoutreix et al. 1695; Ramírez-Salinas et al. 2020). Besides, many available therapeutics to reposition against SARS-CoV-2 infection have been subjected to clinical trials (Hoffmann et al. 2021; Abdool Karim and Oliveira 2021). Nevertheless, the search for other novel molecules cannot pause while waiting for the outcomes because the demand for new therapeutic agents is inflating (Singh et al. 2021).

Since the commencement of the SARS-CoV-2 infection, various strategies have been implemented to tackle the pandemic. One of them consists of drug repositioning using computational techniques (Novac 2013). With this approach, the emergence of commercial medications with therapeutic potential against SARS-CoV-2 infection has been proposed. One such example is the repurposing of imidazole- and benzimidazole-based anti-parasitic drugs. Imidazole and benzimidazole rings are the most important nitrogen-containing heterocycles, which are widely explored and utilized by the pharmaceutical industry for the drug discovery and development. Chemically, imidazole is a five-membered aromatic heterocycle present in the biological building blocks, such as amino acid (histidine), histamine, purines, and biotin (Fei and Zhou 2013), whereas benzimidazole is a six-membered bicyclic heteroaromatic cycle present in vitamin B12 and many other biological compounds (Wright 1951; Barker et al. 1960).

Both the heterocycles contain two nitrogen atoms with amphoteric nature and exist in two equivalent tautomeric forms. In addition, the electron-rich nitrogen heterocycles could readily accept/donate electrons and form diverse weak interactions easily. These unique characteristic features are beneficial for the imidazole and benzimidazole derivatives to readily bind with a variety of therapeutic targets, thereby exhibiting broad range of pharmacological or biological activities (Wright 1951; Bhatnagar et al. 2011; Ingle and Magar 2011; Gaba et al. 2010, 2014; Narasimhan et al. 2011; Yadav and Ganguly 2015). Their ubiquitous properties and key role in management of diseases have attracted special interest in imidazole- and benzimidazole-based medicinal chemistry.

In this paper, we performed docking and binding free energy calculations of available imidazole- and benzimidazole-based anti-parasitic drugs within the catalytic pockets of SARS-CoV-2 Alpha and Gamma variant SGP.

Materials and methods

Preparation of protein

To circumvent the current pandemic of SARS-CoV-2 mutants, X-ray crystallographic protein structures of mutated SGPs, with accession IDs 7NEH (resolution 1.77 Å; R-value work 0.188 Å; R-value free 0.198 Å) (Supasa et al. 2021) and 7NXC (resolution 3.14 Å; R-value work 0.226 Å; R-value free 0.278 Å) (Dejnirattisai et al. 2021), were downloaded from the RCSB protein data bank. Crystal structures were prepared using a protein preparation wizard available in the Maestro interface (Sastry et al. 2013). While preparing, hydrogen atoms were added, and bond orders were assigned using the PROPKA method (Beard et al. 2021). Further, filling of missed side chains and loops was achieved using the module Prime (Iwaloye et al. 2020). Minimization of protein structure was carried out using OPLS3e force field (Roos et al. 2019). A 10 Å grid boxes for the two receptors were generated using the wizard Receptor Grid Generation available in Glide module (Halgren et al. 2004) by selecting the centroid of active amino acid residues Tyr52, Val58, Ser60, Arg61 Gly158, and Asn159 for PDB ID: 7NEH (Supasa et al. xxxx) and Lys26, Asn53, Asn103, Asn343, and Leu368 for PDB ID: 7NXC (Dejnirattisai et al. 2021).

Preparation of ligand

A total of 15 reported imidazole- and benzimidazole-based anti-parasitic drugs were selected for molecular docking studies to discover the potential therapeutics against SARS-CoV-2 mutants. The ligands and standard drugs (hydroxychloroquine and remdesivir) used are given in Table 1. The 3D structures were taken from the PubChem database in .sdf format. Geometry optimization and minimization of the molecules were carried out using the algorithms available in the LigPrep Wizard. Ligands were prepared by adding hydrogens, desalting and ionizing at pH 7 ± 2. OPLS3e force field was used for energy minimization of ligands by keeping root-mean-square deviation (RMSD) cut-off 0.01 Å and standard energy functions of molecular mechanics to produce low-energy conformers (Friesner et al. 2006).

Table 1

Details of the imidazole- and benzimidazole-based anti-parasitic drugs

S. No.	Compound name	PubChem ID	Mol. Wt	Mol. Formula	No. of aHBD	No. of bHBA
1	Albendazole	2082	265.33	C12H15N3O2S	2	4
2	Azanidazole	643,671	246.23	C10H10N6O2	1	6
3	Benznidazole	31,593	260.25	C12H12N4O3	1	4
4	Dimetridazole	3090	141.13	C5H7N3O2	0	3
5	Fenbendazole	3334	299.3	C15H13N3O2S	2	4
6	Flubendazole	35,802	313.28	C16H12FN3O3	2	5
7	Mebendazole	4030	295.29	C16H13N3O3	2	4
8	Megazol	29,698	226.22	C6H6N6O2S	1	7
9	Metronidazole	4173	171.15	C6H9N3O3	1	4
10	Nimorazole	23,009	226.23	C9H14N4O3	0	5
11	Ornidazole	28,061	219.62	C7H10ClN3O3	1	4
12	Pretomanid	456,199	359.26	C14H12F3N3O5	0	9
13	Thiabendazole	5430	201.25	C10H7N3S	1	3
14	Tinidazole	5479	247.27	C8H13N3O4S	0	5
15	Triclabendazole	50,248	359.7	C14H9Cl3N2OS	1	3
16	Hydroxychloroquine	3652	335.9	C18H26ClN3O	2	4
17	Remdesivir	121,304,016	602.6	C27H35N6O8P	4	13

aHydrogen bond donor groups; bHydrogen bond acceptor groups

Molecular docking studies

Molecular docking is an in silico technique used to discover the critical amino acid interactions between the chosen receptors and low-energy conformers prepared using LigPrep. Protein–ligand interactions were categorized by the scoring function, which indirectly demonstrates the ligand’s binding affinity with the chosen receptor. Glide extra-precision (XP) docking protocol was employed without applying any constraints. XP docking was carried out to speculate the ligand efficiency and binding affinity as potent inhibitors of SARS-CoV-2 mutants. Maestro interface was used to visualize the poses of docked complexes.where a and b = coefficient constants of Coulombic and van der Waals energy; Metal = binding with metals; Lipo = lipophilic term; Hbond = hydrogen bonding; RotB = rotatable bonds penalty; BuryP = buried polar groups penalty; Site = catalytic pocket polar interaction.

Binding free energy (BFG) calculations

BFG was calculated to determine the strength of the ligands using molecular mechanics-generalized born surface area (MM-GBSA) technique, Prime module (Jacobson et al. 2004). Energy minimization of protein–ligand docked complex is done using the local optimization feature present in the Prime module. Simulation of the complexes was achieved by input ligand partial charges without smearing any constraints by employing VSGB 2.0 energy model (Li et al. 2011). The energy of the docked complexes was calculated using OPLS3e force field. The BFG (Gbind) was calculated using the following formula:Gcpx = protein–ligand complex free energy; Glig = ligand energy; Grec = unbound receptor energy.

Molecular dynamic (MD) simulations

GROMACS 2018.1 software package was used to perform a 100 ns of MD simulation for flubendazole/7NEH complex. Gromos54a7 force field and PRODRG web server were utilized for the energy optimization and constructing structural topologies of the complex (Schüttelkopf and Aalten 2004). SPC water model and cubic boundary box were employed for the simulation of a system. Neutralization was attained using inorganic ions such as sodium (Na+) and chloride (Cl−) ions, controlling periodic limit constraints. System equilibration was attained at a pressure of 1 atm and temperature of 300 K (Berendsen et al. 1984). Energy minimization of the system was done using the steepest descent algorithm.

The command lines g_gyrate, g_cluster, g_rmsd, and g_rmsf were used for analysing radius of gyration (ROG), RMSD, mean square displacement, and root-mean-square fluctuation (RMSF), respectively. The docked complex flubendazole/7NEH was employed to observe the migration over the simulation time.

Molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) calculation

Binding free energy calculation of flubendazole/7NEH complex was performed using the MM-PBSA approach (Kollman et al. 2000). This technique calculates the binding free energy of the end states directly, by omitting the intermediate states simulation, and integrates the molecular mechanical energies with the continuum solvent strategies. Explicit water molecules were eliminated from the trajectories.

The binding free energy (ΔGbind) was calculated using the following formula:where TΔSconf is the solute entropies calculated by the module NMODE in GROMACS.

ΔEgas indicates the interaction energy between protein and ligand in the gaseous phase:where ΔEelec and ΔEvdW denote the protein–ligand electrostatic and van der Waal interactions, respectively.

Further, total solvation free energy (ΔGsolv) is the sum of a polar (ΔGNP) and electrostatic (ΔGPB) solvation free energies (Luo et al. 2002):

Results and discussion

Molecular docking analysis

To determine a possible drug candidate to combat SARS-CoV-2 mutant strains, molecular docking studies were employed using 15 imidazole- and benzimidazole-based anti-parasitic drugs on the catalytic sites of Alpha and Gamma variant SGP. All these compounds were docked against SARS-CoV-2 Alpha and Gamma variant SPGs and ranked based on their glide score (Gscore). The docking results were compared with standards (hydroxychloroquine and remdesivir) used in the study and are summarized in Tables 2 and 3.

Table 2

Molecular docking (kcal/mol) results of imidazole- and benzimidazole-based anti-parasitic drugs in the catalytic pocket of SARS-CoV-2 Alpha (PDB ID: 7NEH)

Ligands	aGscore	bGvdw	cGcoul	dGenergy	eGmodel	fHBond
Albendazole	− 3.961	− 26.805	− 6.463	− 33.268	− 44.47	− 1.33
Azanidazole	− 3.116	− 24.323	− 4.723	− 29.047	− 35.234	− 1.254
Benznidazole	− 3.03	− 23.111	− 8.65	− 31.762	− 37.55	− 1.33
Dimetridazole	− 2.375	− 12.822	− 3.948	− 16.77	− 20.277	− 1.083
Fenbendazole	− 3.645	− 31.583	− 11.799	− 43.383	− 52.358	− 1.771
Flubendazole	− 4.418	− 24.022	− 14.705	− 38.728	− 55.323	− 1.866
Mebendazole	− 4.066	− 30.258	− 12.084	− 42.342	− 57.592	− 2.463
Megazol	− 2.737	− 22.155	− 6.542	− 28.697	− 32.464	− 0.572
Metronidazole	− 3.216	− 14.119	− 7.66	− 21.779	− 24.783	− 1.866
Nimorazole	− 3.22	− 21.164	− 4.43	− 25.594	− 28.656	− 1.488
Ornidazole	− 3.292	− 17.634	− 7.576	− 25.211	− 31.525	− 1.578
Pretomanid	− 3.197	− 33.478	− 4.254	− 37.732	− 46.419	− 0.834
Thiabendazole	− 2.683	− 21.526	− 3.896	− 25.422	− 31.373	− 0.7
Tinidazole	− 2.367	− 17.938	− 6.987	− 24.926	− 29.456	− 0.567
Triclabendazole	− 3.565	− 25.987	− 4.778	− 30.764	− 39.083	− 0.7
Hydroxychloroquine	− 3.473	− 27.768	− 7.869	− 35.637	− 48.26	− 0.968
Remdesivir	− 6.266	− 41.932	− 13.832	− 55.764	− 67.518	− 3.391

a Glide score; bglide van der Waals energy; cglide Coulombic energy; dglide energy; eglide model energy; fhydrogen bond energy

Table 3

Molecular docking (kcal/mol) results of imidazole- and benzimidazole-based anti-parasitic drugs in the catalytic pocket of SARS-CoV-2 Gamma (PDB ID: 7NXC)

Ligands	aGscore	bGvdw	cGcoul	dGenergy	eGmodel	fHBond
Albendazole	− 1.860	− 21.685	− 2.822	− 24.508	− 27.891	− 0.002
Azanidazole	− 1.594	− 21.520	− 1.994	− 23.514	− 28.139	− 0.700
Benznidazole	− 2.373	− 26.945	− 3.686	− 30.631	− 37.892	− 0.700
Dimetridazole	− 1.829	− 16.753	− 2.451	− 19.205	− 21.913	0.000
Fenbendazole	− 3.148	− 21.431	− 8.278	− 29.709	− 37.828	− 0.748
Flubendazole	− 2.258	− 21.343	− 6.254	− 27.597	− 31.552	− 0.574
Mebendazole	− 4.162	− 24.371	− 5.172	− 29.542	− 37.389	− 1.243
Megazol	− 2.458	− 18.779	− 8.69	− 27.468	− 31.91	− 0.35
Metronidazole	− 3.214	− 25.075	− 1.051	− 26.126	− 32.506	0
Nimorazole	− 1.68	− 19.485	− 4.419	− 23.904	− 31.013	− 0.35
Ornidazole	− 2.933	− 17.969	− 7.322	− 25.291	− 28.81	− 0.7
Pretomanid	− 3.3	− 25.233	− 4.805	− 30.038	− 40.936	− 0.459
Thiabendazole	− 2.562	− 20.363	− 1.696	− 22.06	− 23.979	0
Tinidazole	− 1.251	− 11.736	− 6.678	− 18.414	− 20.838	− 0.191
Triclabendazole	− 1.567	− 22.102	− 5.012	− 27.114	− 32.911	− 0.700
Hydroxychloroquine	− 2.536	− 22.307	− 5.582	− 27.889	− 33.61	− 0.591
Remdesivir	− 3.115	− 34.103	− 5.997	− 40.1	− 47.39	− 0.7

a Glide score; bglide van der Waals energy; cglide Coulombic energy; dglide energy; eglide model energy; fhydrogen bond energy

Interaction analysis of ligands with SARS-CoV-2 Alpha variant SGP

All the compounds exhibited significant activity in the range of − 2.367 to − 4.418 kcal/mol in the catalytic pocket of SARS-CoV-2 Alpha variant SGP (Table 2 and Supplementary Figure S1a-q). Gratifyingly, all the imidazole- and benzimidazole-based anti-parasitic drugs showed excellent results in computational docking studies. Glide score exhibited a good correlation with glide energy (R2 = 0.8541), glide emodel (R2 = 0.8482), and Hbond (R2 = 0.8398), while the moderate correlation with glide van der Waal (R2 = 0.7601) and glide Coulomb (R2 = 0.7284). Further, the ligands exhibited significant interactions at the same catalytic site compared to the standard drugs hydroxychloroquine and remdesivir. Hydroxychloroquine exhibited significant interactions with the Ser349 and Tyr351 with a Gscore of − 3.471 kcal/mol (Supplementary Figure S1p), whereas remdesivir has exhibited significant interactions with amino acid residues Ser349, Ala352, Asn354, and Asn450 with a Gscore of − 6.266 kcal/mol (Supplementary Figure S1q).

Out of 15 test compounds, flubendazole exhibited more potent inhibitory activity as indicated with its high Gscore of − 4.481 kcal/mol. The other compound, mebendazole, also showed a promising Gscore of − 4.066 kcal/mol within the catalytic pocket of protein structure.

Flubendazole has shown five hydrogen bonding (HB) interactions similar to that of remdesivir with the critical interacting residues of the catalytic pocket of SARS-CoV-2 Alpha lineage SGP (Fig. 1). The oxygen of the terminal ester attached at position two of benzimidazole formed a HB network with the protonated amino group of Asn354 ((> C = O)O⋯HN(H), 2.41 Å), whereas -OH and -NH of Ser349 residue involved in the formation of two HB interactions with the carbonyl group of carbamate linker attached to benzimidazole group at position two (> C = O⋯HN(H), 1.85 Å; > C = O⋯HO, 2.09 Å). The carbonyl group between the two aromatic rings has shown two HB interactions with the -NH moieties of Asn448 (> C = O⋯HN(H), 1.79 Å) and Asn450 (> C = O⋯HN(H), 2.00 Å), respectively. Moreover, flubendazole is stabilized by two π–π stacking interactions, observed between electron clouds of the benzimidazole nucleus and the phenyl ring of Tyr451.

Fig. 1

3D interaction diagram of flubendazole with amino acid residues necessary for the inhibition of SARS-CoV-2 Alpha variant

Further, the hydrophilic–hydrophobic map of the flubendazole reveals that the nitrogen moieties present in the benzimidazole ring and HNO = C(O)CH3 fragment are located well in the hydrophilic region. In contrast, phenyl rings are buried in the hydrophobic part of the receptor. On the other hand, the fluorobenzene group at position five of the benzimidazole ring is well buried in the hydrophilic pocket because of the high dielectric constant of fluorinated compounds compared to the simple aromatic rings (Biffinger et al. 2004) (Supplementary Figure S2).

Interaction analysis of ligands with SARS-CoV-2 Gamma variant SGP

Comparison of the docked poses of conformationally flexible ligands along with the standards revealed similar kinds of binding poses within the catalytic pocket of SARS-CoV-2 Gamma variant SGP (Table 3 and Supplementary Figure S3). Glide score of selected analogues in the catalytic pocket of SARS-CoV-2 Gamma variant SGP was observed in the range of − 1.251 to − 4.162 kcal/mol. The critical amino acid interactions were dominated in Phe347-Ser357, Tyr396-Pro426, Asn450-Leu475, and Ser514-Glu516. Glide score exhibited a good correlation with glide emodel (R2 = 0.6077), while the moderate correlation with glide van der Waal (R2 = 0.5853) and glide energy (R2 = 0.5631). The results obtained were compared with that of standards; hydroxychloroquine exhibited significant interactions with the Thr430 and Glu516 with a Gscore of − 2.536 kcal/mol (Supplementary Figure S3p). In comparison, remdesivir has exhibited significant interactions with amino acid residues Arg355, Phe464, Arg466, and Glu516 with a Gscore of −3.115 kcal/mol (Supplementary Figure S3q) (Table 3).

Table 4

Binding free energy (kcal/mol) of imidazole- and benzimidazole-based anti-parasitic drugs in the binding pocket of SARS-CoV-2 Alpha variant SGP (PDB ID: 7NEH) enzyme calculated using MM-GBSA approach

Ligand	aΔGbind	bΔGcoul	cΔGcov	dΔGHbond	eΔGlipo	fΔGsolvGB	gΔGvdw
Albendazole	− 50.01	− 28.09	2.29	− 4.56	− 10.61	22.87	− 28.3
Azanidazole	− 26.81	12.49	3.26	− 1.36	− 6.43	− 0.06	− 29.39
Benznidazole	− 43.16	− 15.6	− 6.84	0.09	− 4.69	5.87	− 17.12
Dimetridazole	− 24.31	− 8	− 5.63	0.38	− 2.6	5.05	− 13.57
Fenbendazole	− 57.54	− 48.46	0.55	− 1.55	− 8.82	33.55	− 27.15
Flubendazole	− 66.05	− 40.26	6.41	− 6.83	− 12.21	20.92	− 45.62
Mebendazole	− 55.22	− 17.88	2.22	− 1.25	− 9.82	21.74	− 43.03
Megazol	− 36.55	− 2.41	− 4.56	− 3.67	− 5.81	5.2	− 21.06
Metronidazole	− 36.1	− 22.15	2.52	− 2.93	− 6.14	4.33	− 10.15
Nimorazole	− 36.47	− 0.89	− 3.49	0.09	− 8.84	− 4.02	− 17.54
Ornidazole	− 35.21	− 30.32	6.66	− 4.25	− 6.1	15.48	− 15.77
Pretomanid	− 63.77	− 2.43	− 4.54	− 4.12	− 15.18	2.95	− 39.05
Thiabendazole	− 45.2	− 24.03	− 1.39	− 3.6	− 8.75	19.41	− 19.73
Tinidazole	− 13.63	36.06	− 3.56	− 1.58	− 6.38	− 13.09	− 24.93
Triclabendazole	− 38.4	− 4.97	− 12.96	1.3	− 13.26	25.48	− 32.24
Hydroxychloroquine	− 59.25	− 23.76	− 3.8	− 2.67	− 13.66	13.5	− 27.68
Remdesivir	− 49.95	− 18.81	− 0.73	− 2.05	− 13.48	31.11	− 44.89

aFree energy of binding; bCoulombic energy; cCovalent energy (internal energy); dhydrogen bonding energy; ehydrophobic energy (non-polar contribution estimated by solvent accessible surface area); felectrostatic solvation energy; gvan der Waals energy

Amongst 15 test ligands, mebendazole exhibited more potent inhibitory activity as indicated with its high Gscore of − 4.162 kcal/mol (Fig. 2). The other compounds showed moderate biological activity within the binding pocket of protein structure. Mebendazole has shown three HB interactions with the critical interacting residues of SARS-CoV-2 Gamma variant SGP. The NH moiety of the benzimidazole nucleus exhibited a HB interaction with the carbonyl group of Phe347 (NH⋯O = C < , 1.92 Å). The Ser349 formed HB two interactions with the carbonyl group of carbamate fragment (> C = O⋯HN, 1.85 and 2.05 Å). Further, the carbamate fragment’s NH moiety showed a HB interaction with the carbonyl group of Asn450 (NH⋯O = C < , 2.57 Å).

Fig. 2

3D interaction diagram of mebendazole with amino acid residues necessary for the inhibition of SARS-CoV-2 Gamma variant

In addition, the hydrophobic–hydrophilic map of the mebendazole demonstrated that the imidazole ring of the benzimidazole scaffold is well located in the hydrophilic region of the catalytic pocket, while the aromatic group is present at position five of the benzimidazole ring. At the same time, the carbonyl group present between the two aromatic rings is well buried in the hydrophilic pocket of the receptor (Supplementary Figure S4).

Binding free energy calculation analysis

In the current study, the BFG of selected 15 molecules was calculated using the MM-GBSA strategy to get insight into the relative potencies and binding strengths against SARS-CoV-2 Alpha and Gamma variant SGP.

Binding affinity and free energy analysis of ligands against SARS-CoV-2 Alpha variant SGP

The predicted BFG score of docked poses of ligands was found to be in the range of − 24.31 to − 66.05 kcal/mol (Table 4). Similar results were observed for the standard compounds hydroxychloroquine (ΔGbind = -59.25 kcal/mol) and remdesivir (ΔGbind = -49.95 kcal/mol). It is evident from the study that van der Waal (ΔGvdW = − 10.15 to – 43.03 kcal/mol) energy term is the major contributor to inhibitory activity in the catalytic pocket of 7NEH, while Coulombic energy term (ΔGcou = 36.06 to – 48.46 kcal/mol) moderately favours the ligand binding in most of the active ligands. In all the active compounds, non-polar solvation (ΔGlipo = − 2.6 to – 15.18 kcal/mol) and covalent energy term are considered unfavourable for the ligand’s inhibitory activity within the binding site of 7NEH. Nevertheless, the compound flubendazole showed relatively BFG (ΔGbind = − 66.05 kcal/mol), van der Waal energy (ΔGvdW = − 45.62 kcal/mol), and Coulombic energy (ΔGcou = − 40.26 kcal/mol). It is evident from the result that non-polar van der Waals (ΔGvdW) is considered to be the main contributor to the inhibitory activity in the binding pocket of 7NEH.

Binding affinity and free energy analysis of ligands against SARS-CoV-2 Gamma variant SGP

The results of BFG of selected compounds against SARS-CoV-2 Gamma variant SGP are represented in Table 5. Similar results were observed for the standard compound drugs hydroxychloroquine (ΔGbind = − 66.43 kcal/mol) and remdesivir (ΔGbind = − 75.45 kcal/mol). The predicted BFG was observed to be in the range of − 22.35 to − 71.54 kcal/mol. It is clear from the studies that van der Waals (ΔGvdW = − 10.15 to − 43.03 kcal/mol) and Coulombic energy (ΔGcou = − 10.15 to − 43.03 kcal/mol) terms are found to be highly favourable for the inhibitory activity of the ligands in the active site of 7NXC, whereas non-polar solvation energy terms and covalent energy terms fluctuate between 19.66 to − 5.83 kcal/mol and − 12.08 to − 6.16 kcal/mol, clearly showing the low preference of their components in the ligand binding. Nonetheless, the compound mebendazole showed relatively high BFG (ΔGbind = − 71.54 kcal/mol), van der Waal energy (ΔGvdW = − 29.7 kcal/mol) and Coulombic energy (ΔGcou = − 25.51 kcal/mol). Thus, it is evident that van der Waals term and Coulombic energy terms are the key driving forces for the inhibitory activity of these ligands.

Table 5

Binding free energy (kcal/mol) of imidazole- and benzimidazole-based anti-parasitic drugs in the binding pocket of SARS-CoV-2 Gamma variant SGP (PDB ID: 7NXC) enzyme calculated using MM-GBSA approach

Ligand	aΔGbind	bΔGcoul	cΔGcov	dΔGHbond	eΔGlipo	fΔGsolvGB	gΔGvdw
Albendazole	− 46.97	− 37.53	2.02	− 3.04	− 2.82	17.14	− 15.71
Azanidazole	− 28.84	− 4.03	4.41	− 0.96	− 3.54	2.40	− 23.16
Benznidazole	− 50.75	− 25.62	6.92	− 3.36	− 11.34	19.66	− 32.5
Dimetridazole	− 36.08	− 9.52	− 9.79	− 0.43	3.77	− 7.04	− 18.42
Fenbendazole	− 38.13	− 33.76	3.16	− 2.59	− 11.46	27.47	− 16.12
Flubendazole	− 48.84	− 38.1	10.97	− 2.38	− 10.59	19.8	− 23.29
Mebendazole	− 71.54	− 25.51	− 6.16	− 3.74	− 11.81	7.72	− 29.7
Megazol	− 37.02	9.79	− 3.9	− 1.11	− 5.7	− 5.66	− 27.37
Metronidazole	− 18.82	13.88	6.62	1.85	− 7.75	− 5.83	− 24.65
Nimorazole	− 36.12	− 1.85	2.17	− 3.31	− 12.3	0.64	− 22.58
Ornidazole	− 28.6	− 6.99	10	0.56	− 2.53	− 3.01	− 23.4
Pretomanid	− 43.1	− 17.85	8.19	− 1.42	− 11.08	11.77	− 26.68
Thiabendazole	− 43.03	− 32.14	12.08	− 2.33	− 8.5	13.06	− 16.88
Tinidazole	− 22.35	1.3	2.61	− 0.58	− 1.17	− 5	− 19.55
Triclabendazole	− 26.02	− 24.18	− 0.61	− 0.00	− 5.39	19.15	− 15.09
Hydroxychloroquine	− 66.43	− 25.53	10.38	− 5.99	− 19.58	7.27	− 30.84
Remdesivir	− 75.45	− 8.01	− 0.35	− 3.25	− 27.42	12.49	− 42.27

aFree energy of binding; bCoulombic energy; cCovalent energy (internal energy); dhydrogen bonding energy; ehydrophobic energy (non-polar contribution estimated by solvent accessible surface area); felectrostatic solvation energy; gvan der Waals energy

MD simulation analysis

To investigate the docking validity, the docked complex of flubendazole/7NEH was subjected to 100 ns MD simulation. RMSD provides data related to the protein’s backbone, and the value was found to be 3.6 nm; small initial fluctuation was observed, which was stabilized after 10 ns due to the decrease in amplitude of fluctuation with time (Fig. 3).

Fig. 3

RMSD plot of backbone atoms of flubendazole/7NEH over a time period of 100 ns MD simulation

In addition, RMSF was computed to determine the dynamic behaviour of atomic positional fluctuations of backbone residues and the value was found to be 0.45 nm at 200–300 residues (Fig. 4).

Fig. 4

RMSF plot of backbone atoms of 7NEH over a time period of 100 ns MD simulation

ROG was computed to determine the compactness of the complex and protein structure, and the value was observed between the range of 1.42 to 1.49 nm (Fig. 5).

Fig. 5

Radius of gyration (ROG) plot of flubendazole/7NEH complex for 100 ns MD simulation

MM-PBSA

MM-PBSA analysis aids in segregating the total binding free energy into electrostatic, van der Waals, and solute–solvent interactions, to get insight into the binding modes and protein–ligand interaction process. The total binding free energy of the complex was found to be -68.831 kcal/mol. From the results, it is evident that van der Waals interactions (−79.735 kcal/mol) play an essential role in the MD simulation, contributing noticeably more to total interaction energy than electrostatic interactions (-55.919 kcal/mol) indicating dominating effect of polar components in the system. The polar/electrostatic and apolar/hydrophobic contributions to free energy are also recorded in Table 6. Further, SASA was calculated to determine the hydrophobic core in the protein structure and the value lies between 80 and 95 nm2, indicating the stability throughout the MD simulation (Fig. 6). SASA has not exhibited any fluctuations after attaining equilibrium throughout the MD simulation study. In disparity, the ΔEPB and SASA energy display limited contribution to the energy.

Table 6

MM-PBSA analysis of compound flubendazole/7NEH complex

ΔEvdW	ΔEelec	ΔEPB	SASA energy	Binding energy
− 79.735 ± 8.449 kcal/mol	− 55.919 ± 19.029 kcal/mol	104.009 ± 19.002 kcal/mol	− 59.004 ± 0.789 kcal/mol	-68.831 ± 19.327 kcal/mol

Fig. 6

Solvent accessible surface area (SASA) of flubendazole/7NEH complex

Conclusion

Several studies are ongoing using small molecule inhibitors to circumvent the life-threatening infections caused by the mutant strains of SARS-CoV-2. Imidazole- and benzimidazole-based anti-parasitic drugs have been suggested as effective treatments in tackling SARS-CoV-2 mutants. In this study, molecular docking studies were performed on 15 imidazole- and benzimidazole-based anti-parasitic drug molecules along with the standard drugs hydroxychloroquine and remdesivir. Amongst these compounds, flubendazole and mebendazole are useful as potential drug candidates against SARS-CoV-2 mutants, as evidenced by the glide score (− 4.418 kcal/mol for 7NEH and − 4.162 kcal/mol for 7NXC kcal/mol) and BFG (− 66.05 kcal/mol for 7NEH and − 71.54 kcal/mol for 7NXC). Further, flubendazole/7NEH complex was subjected to MD simulation and the results obtained satisfied all the parameters of trajectory analysis. We, therefore, propose that the inhibitory potentials of flubendazole and mebendazole against spike glycoprotein variants of SARS-CoV-2 should further be explored in the subsequent wet laboratory experiments. This study facilitates the discovery process of potent inhibitors against SARS-CoV-2 mutants.

Below is the link to the electronic supplementary material.

Supplementary file1 (PDF 4636 kb)