Repurposing benzimidazole and benzothiazole derivatives as potential inhibitors of SARS-CoV-2: DFT, QSAR, molecular docking, molecular dynamics simulation, and in-silico pharmacokinetic and toxicity studies

1 Introduction

The current COVID-19 pandemic, caused by the highly contagious newly emerged SARS-CoV-2, is a global health threat that affects millions of people worldwide, and impacts negatively on all sectors of society (Mohapatra et al., 2020a; b; Mohapatra and Rahman, 2021). Although the fatality rate for SARS-CoV-2 (∼4%) is lower compared to that of SARS-CoV (∼10%) and MERS-CoV (∼35%), its transmission rate is much higher than that of previously reported coronaviruses (CoVs) (Guo et al., 2020). In many countries, this led to lockdown measures being put into place to minimize human-to-human transmission (Mohapatra et al., 2021b). COVID-19 related deaths have now surpassed those from other epidemics caused by pathogens such as influenza A & B, Zika, Ebola and other CoVs (Mahal et al., 2021) and declared global pandemic (Dhama et al., 2020a).

The SARS-CoV-2 virus is mainly transmitted through respiratory aerosols/droplets, and to some extent via the fecal-oral route (Chan et al., 2020; Sah et al., 2021). Other means of transmission include direct and indirect contacts with contaminated objects/surfaces/fomites (Mohapatra et al., 2021a; Morawska and Cao, 2020). The SARS-CoV-2 virus primarily targets the upper and lower respiratory tract, and most people infected will experience only mild or moderate symptoms. However, in some individuals such as the elderly and/or those with pre-existing comorbidities (e.g. obesity, diabetes, cardiovascular disorders), SARS-CoV-2 infection can lead to severe complications, with multiple organ damage (Dhama et al., 2020b; Arumugam et al., 2020).

Scientists throughout the world continue to endeavor to discover and develop suitable and effective approaches to tackle this pandemic. This includes the development of vaccines such as the Covishield, Covaxin (India), AstraZeneca (UK), Moderna, Novavax (USA), Pfizer/BioNTech (Germany), Sputnik V (Russia), CanSino, Sinovac and BBIBP-CorV (China) vaccines (Madkaikar et al., 2021; Bagcchi, 2021). Vaccination programs are underway in many countries; even if clinical efficacy against the new SARS-CoV-2 variants has yet to be established for all vaccines (Callaway, 2021). Efforts of the scientific community have also been focused on the discovery of new antiviral drug leads. Several reports have highlighted the importance of targeting the M^pro of SARS-CoV-2 - an enzyme with a pivotal role in mediating viral replication and transcription - and have used molecular docking as a promising strategy to design and develop possible M^pro inhibitors (Ancy et al., 2020; Baildya et al., 2020; Mohapatra et al., 2021c). Another aim has been to target virus entry mechanisms via disruption of the interaction of SARS-CoV-2 with the human ACE2 receptor responsible for viral attachment to the host cells (Rezaei et al., 2021; Joshi et al., 2020).

In this work, we used a comprehensive in-silico approach, including molecular docking and MD simulation studies, to investigate the predictive binding affinity of four benzimidazoles and two benzothiazoles previously synthesized by our group (Mohapatra et al., 2011; 2018) towards SARS-CoV-2 M^pro and the human ACE2 receptor. We also report on the pharmacokinetic, drug-likeness, and toxicological properties of the studied compounds.

2 Density functional theory (DFT) studies

DFT is a computational modeling study used to investigate the electronic structure of molecules. The DFT analysis was carried out through GAUSSIAN inter phase. All compounds were first optimized using the Gauss View 6.0.16 program (Gauss View 6.0, 2019), and then subjected to DFT calculation using the GAUSSIAN 09 suite programs (Frisch et al., 2009). The computation was performed using the B3LYP method and exchange–correlation functional with a 6–31 + G (d, p) basis set for carbon, nitrogen, oxygen, sulfur, and hydrogen atoms. Relevant energetic properties such as the dipole moment (D) and the single point energy were calculated for each compound. Comparison of the calculated energy values revealed that compound 4 had greater single point energy compared with the other compounds (El‐ajaily et al., 2019; Sarangi et al., 2018; 2020; Mahapatra et al., 2013). Compound 2 had a lower energy value which indicated larger stability, while compound 3 possessed a greater D value compared to other compounds (Table S1).

2.1

2.1 Frontier molecular orbitals (FMO) studies

FMO studies can predict the chemical reactivity of compounds and identify their most likely reactive sites. The calculated energies E_HOMO and E_LUMO of compounds help to explain global reactivity descriptors such as chemical hardness, chemical potential, and electrophilicity. The stability of the investigated compounds was confirmed from the negative values obtained for their E_HOMO and E_LUMO (Yousef et al., 2013). The electron cloud of compound 1 was localized on the phenolic (–OH) group and azomethine nitrogen (C = N) group relating to the participating atomic sites involved in nucleophilic attack. The band energy gap correlates with the chemical reactivity and chemical stability of molecules. The energy difference [E_HOMO - E_LUMO] for compound 3 was found to be less than the band energy gap observed for other compounds, establishing its greater reactivity. Compound 1, on the other hand, had a greater stability. The band energy values and E_HOMO - E_LUMO energy variation suggest charge transfer (CT) interaction. In this work, the chemical reactivity parameters such as electronegativity (χ), chemical potential (μ), global electrophilicity index (ω), global hardness (η) and global softness (S) were calculated (Table S2). We observed that the chemical softness (S) of compound 4 was lesser than the other compounds, indicating a higher stability for compound 4. Another important parameter is the electrophilicity ( $\omega$ ) of compounds, which assesses its ability to accept electron(s) from its surrounding. Compound 1 showed less electrophilicity than other compounds (Table S2) (Fig. 1). $$\chi =-\frac{({\mathbf{E}}_{\mathbf{L}\mathbf{U}\mathbf{M}\mathbf{O}}+{\mathbf{E}}_{\mathbf{H}\mathbf{O}\mathbf{M}\mathbf{O}})}{2}$$ $$\mu =-\chi =\frac{({\mathbf{E}}_{\mathbf{L}\mathbf{U}\mathbf{M}\mathbf{O}}+{\mathbf{E}}_{\mathbf{H}\mathbf{O}\mathbf{M}\mathbf{O}})}{2}$$ $$\eta =\frac{({\mathbf{E}}_{\mathbf{L}\mathbf{U}\mathbf{M}\mathbf{O}}-{\mathbf{E}}_{\mathbf{H}\mathbf{O}\mathbf{M}\mathbf{O}})}{2}$$ $$\mathbf{S}=\frac{1}{2\eta }$$ $$\omega =\frac{{\mu }^2}{2\eta }$$ $$\sigma =\frac{1}{\eta }$$

Close

Fig. 1

Optimized structures of compounds 1–6 using FMOs.

Export to PPT

3 QSAR studies

QSAR modelling was used to characterize the physico-chemical properties of the investigated compounds and predict their correlations to biological activity. QSAR studies were performed using the Hyper Chem Professional 8.0.3 program. First, the structure of each compound was optimized with a semi-empirical PM3 method using a (MM⁺) force field, and then energy minimization was achieved using a Fletcher-Reeves conjugate gradient algorithm. A range of QSAR parameters including free energy, hydration energy, total energy, refractivity, surface area, volume, polarizability, mass, dipole moment, partition coefficient (log P), and RMS gradient were calculated. Interestingly, the log P value for compound 5 was the highest amongst all compounds, indicating the ability for this compound to efficiently permeate biological membranes (Mohapatra et al., 2021c) (Table 1).

4 Molecular docking study

The docking analysis was executed using a receptor-oriented flexible docking protocol with the Autodock 4.2.6 software package and Python 3.8.2. The ligands were generated with MGL Tools 1.5.6 and optimized using the Avogadro program. The 3D X-ray structure for SARS-CoV-2 M^pro and the ACE2 receptor was retrieved from the Protein Data Bank ( http://www.rcsb.org/pdb) and used for the molecular docking (Adem et al., 2020; Shivanika et al., 2020; Khan et al., 2020). The active site of each enzyme was fixed with Meta Pocket 2.0 and the docking cavities were selected for the target proteins. All ligand interactions with the protein amino acid residues were visualized using the Discovery Studio 3.5 software.

The results of docking analyses were expressed in terms of scoring function (free binding energy, ΔG) for the best conformational positions of each ligand for comparative analysis. The binding affinity of the investigated compounds against SARS-CoV-2 M^pro and ACE2 is related to the favorable geometric arrangement of these ligands in the active site, hydrogen bond formation, donor–acceptor interactions and intermolecular electrostatic interactions which ultimately leads to the formation of stable complexes. The thermodynamic probability of the binding was established for each compound by the negative binding energy (ΔG) values. We observed that compounds (1,2,3) showed best results against SARS-CoV-2 M^pro and compounds (1,2,3,5) showed best results against the human ACE2 receptor (Tables S3 & S4).

The active site of M^pro (PDBID: 6LU7) is occupied by charged, polar and hydrophobic amino acid residues. A detailed description of the ligand–protein interactions within the active site of M^pro is presented in Table S3. Compound-1 formed H-bonds with Thr199 (OG1), Leu287 (O), and ionic interactions with Asp289 (OD1, OD2) (Fig. 2). Interactions with Thr199, Leu287 and Asp289 were observed for all compounds. When docked with the human ACE2 receptor (PDB ID: 6 M18), compound-1 showed ionic interactions with Asp494 (OD1), H-bonding with His493, and pi-pi interaction with Arg671 (Table S4). Interactions with His493, Asp494 and Arg671 were observed for almost all compounds.Fig. 3.

Close

Fig. 2

Docked pose of compounds (1–6) in the SARS-CoV-2 M^pro binding site, and corresponding 2D plots of molecular interactions.

Export to PPT

Close

Fig. 3

Docked pose of compounds (1–6) in the ACE2 binding site, and corresponding 2D plots of molecular interactions.

Export to PPT

5 Molecular dynamics simulation

5.1

5.1 Protein RMSD

MD simulations were performed to study the physical movements and dynamic evolution of the (Comp.1-6LU7) and (Comp.1–6 M18) complexes (Adem et al., 2020; Shivanika et al., 2020). The MD simulation was run using the Schrödinger’s Desmond package following a previously reported procedure (Pant et al., 2020). Ligand-protein complexes were simulated for 100 ns to analyze the RMSD values, which estimate the complex stability (with lower RMSD values indicating greater stability). The protein RMSD trajectory for SARS-CoV-2 M^pro and ACE2 fluctuated up to 5 and 15 ns, respectively, before gradually reaching an equilibrium. The (Comp.1-6LU7) complex completed the 100 ns simulation without much fluctuation, suggesting a greater stability. The RMSD average value of M^pro (6LU7) and ACE2 (6 M18) at equilibrium was found to be 1.8 and 9.5 Å, respectively (Fig. 4). The deviation within 1–3 Å is acceptable for small globular protein. Based on these results, the (Comp.1-6LU7) and (Comp.1–6 M18) complexes were further analyzed for protein–ligand contacts.

Close

Fig. 4

Protein RMSD trajectory of (A) (Comp.1-6LU7) and (B) (Comp.1-6 M18) complex.

Export to PPT

5.2

5.2 Ligand properties

The ligand properties were evaluated by calculating Ligand RMSD, Molecular Surface Area (MolSA), Radius of Gyration (rGyr), Polar Surface Area (PSA) and Solvent Accessible Surface Area (SASA) values. The root-mean-square deviation (RMSD) values fluctuated until 5 and 10 ns for Comp.1-6LU7 and Comp.1–6 M18, respectively, before gradually reaching an equilibrium. The RMSD values were around 0.8 to 2.1 Å (equilibrium ∼ 1.7 Å), and 1.1 to 1.9 Å (equilibrium ∼ 1.3 Å) for Comp.1-6LU7 and Comp.1–6 M18, respectively. The rGyr values fluctuated slightly until 5 and 10 ns for Comp.1-6LU7 and Comp.1–6 M18, respectively, before gradually reaching equilibrium. The rGyr values were around 3.57 to 4.10 Å (equilibrium ∼ 3.90 Å), and 3.53 to 4.35 Å (equilibrium ∼ 3.87 Å) for Comp.1-6LU7 and Comp.1–6 M18, respectively. The MolSA values (equivalent to van der Waals surface area) were constant throughout the simulation. These were around 290 to 300 Å² (equilibrium ∼ 295 Å²), and around 285 to 295 Å² (equilibrium ∼ 293 Å²) for Comp.1-6LU7 and Comp.1–6 M18, respectively. The SASA values fluctuated slightly until 5 and 25 ns for Comp.1-6LU7 and Comp.1–6 M18, respectively, before gradually reaching an equilibrium. These were around 160 to 320 Å² (equilibrium ∼ 210 Å²), and around 60 to 240 Å² (equilibrium ∼ 80 Å²) for Comp.1-6LU7 and Comp.1–6 M18, respectively. The PSA values were constant throughout the simulation. These were around 84 to 102 Å² (equilibrium ∼ 96 Å²), and around 88 to 112 Å² (equilibrium ∼ 96 Å²) for Comp.1-6LU7 and Comp.1–6 M18, respectively (Figs. 5 & 6). In both cases, the ligand showed minimal fluctuation in the initial part of the simulation period, then gradually reached equilibrium which remained constant throughout the simulation, suggesting they possessed good stability in the active site of both proteins.

Close

Fig. 5

Ligand property trajectory for the (Comp.1-6LU7) complex.

Export to PPT

Close

Fig. 6

Ligand property trajectory for the (Comp.1–6 M18) complex.

Export to PPT

5.3

5.3 Protein-ligand contacts study

The protein–ligand (P-L) contacts of the (Comp.1-6LU7) and (Comp.1–6 M18) complexes were studied using contact histograms (Dariya and Nagaraju, 2020). The P-L interactions were classified into four types such as hydrophobic, ionic, water bridges, and hydrogen bonds. The latter are commonly grouped into four sub-types such as side-chain acceptor, side-chain donor, backbone donor and backbone acceptor. They play a crucial role in ligand–protein binding and are essential to take into account when it comes to the design of new drugs. The histograms identified Tyr-239, Leu-271, Met-276, Asn-277, Leu-287, and Asp-471, Lys-475, Glu-489, Pro-492, His-493, Asp-494, Asp-637, Val-672, Asn-674, Leu-675 as the H-bond interacting amino acid residues for the (Comp.1-6LU7) and (Comp.1–6 M18) complexes, respectively. The hydrophobic contacts, involving a hydrophobic amino acid and an aromatic or aliphatic group on the ligand, were characterized as pi-cation, pi-pi and other non-specific interactions. The amino acids that showed hydrophobic interactions with the ligand included Tyr-237, Leu-272, Met-276, Leu-286, Leu-287 for 6LU7, and Trp-478, Pro-492, Met-640, Val-672, Leu-675 for 6 M18. We observed that Asp-494 showed minimal ionic interactions with the ligand for 6 M18. Most of the interacting amino acids (Tyr-239, Met-276, Asn-277, Ala-285, Leu-287 for 6LU7 and Asp-494, Leu-675, Lys-475, Val-672, Asn-674 for 6 M18) showed water bridges with the ligand (Figs. 7 & 8). An overview of all protein–ligand interactions is summarized in Fig. 7(c) & 8(c) - with the top panel showing the total number of specific protein–ligand contacts formed, and the bottom panel showing the nature of the interacting amino acid residues over the course of the trajectory. The residues that made a number of specific contact with the ligand in a particular trajectory frame are represented by a darker shade of orange. The number of P-L contacts varied between 0 and 6 for 6LU7 and between 0 and 12 for 6 M18. The receptor-ligand complex showed darker bands for the Leu-287, Met-276, Tyr-239 row for Comp.1-6LU7 and the Leu-675 and Asp-494 row for Comp.1–6 M18, suggesting that these amino acids have more interactions with the ligand in all possible orientations, in agreement with the histogram results.

Close

Fig. 7

Contact plots and ligand–protein interaction residues for the (Comp.1-6LU7) complex.

Export to PPT

Close

Fig. 8

Contact plots and ligand–protein interaction residues for the (Comp.1-6 M18) complex.

Export to PPT

6 Pharmacokinetics, drug-likeness, and toxicological properties

Knowledge of the pharmacokinetic and drug-likeness properties of newly designed ligands for specific targets is essential to predict the absorption, distribution, metabolism, and excretion of potential new drug leads (Shivanika et al., 2020; Savale et al., 2021). In this study, pharmacokinetic properties predictions for the investigated compounds (Fig. 9) were carried out using the admetSAR (versions 1 and 2) and SwissADME software. Their Lipinski’s properties were determined using physico-chemical data retrieved from the PubChem database. We observed that all compounds obeyed the Lipinski’s rule of five, with molecular weights < 500 g/mol, number of hydrogen bond donors and acceptors < 5, logP value < 5 and molar refractivity < 140 (Lipinski, 2000, 2004). The Topological Polar Surface Area (TPSA) of all compounds was < 90 Å², indicating their ability to permeate biological membranes (Shivanika et al., 2020; Dariya and Nagaraju, 2020). Moreover, all compounds showed good solubility and high gastrointestinal absorption. Their low number of rotatable bonds (3 or 4) suggested a good flexibility. The synthetic accessibility value of all compounds was ≤ 3, which indicated the feasibility of synthesis (Tables S5-S11). The toxicity of the investigated compounds was also evaluated using the pkCSM software tool (http://structure.bioc.cam.ac.uk/pkcsm) (Table S10). We observed that the compounds were neither hepatotoxic nor carcinogenic, and that compound 4 and 6 showed better toxicity results.

Close

Fig. 9

Pharmacokinetic and drug-likeness properties of the investigated compounds (LIPO: lipophilicity, INSOLU: insolubility, INSATU: insaturation, and FLEX: flexibility).

Export to PPT

7 Conclusion

Four benzimidazole and two benzothiazole derivatives were investigated for their binding affinity for M^pro and ACE2 receptor using a molecular docking approach. Compounds 1, 2, and 3 were found to bind with equal affinity to both targets, which may be more effective in clinical practice. One of the benzimidazoles (compound 1) was further subjected to molecular dynamics simulation to explain protein–ligand complex stability, ligand properties, and protein–ligand contacts. We established that all compounds obeyed the Lipinski’s rule of five and showed good solubility, high gastrointestinal absorption, good ability to permeate biological membranes, good synthetic accessibility, and low toxicity. Taken altogether, these data suggest that the investigated compounds have the potential to combat SARS-CoV-2. We hope that this work will provide scope for further studies to establish the potential of benzimidazole and benzothiazole derivatives as new M^pro and/or ACE2 receptor inhibitors.