key: cord-0980321-lpjtdmij
authors: Srivastava, Ritika; Gupta, Sunil K.; Naaz, Farha; Sen Gupta, Parth Sarthi; Yadav, Madhu; Singh, Vishal Kumar; Singh, Anuradha; Rana, Malay Kumar; Gupta, Satish Kumar; Schols, Dominique; Singh, Ramendra K.
title: Alkylated benzimidazoles: Design, synthesis, docking, DFT analysis, ADMET property, molecular dynamics and activity against HIV and YFV
date: 2020-10-06
journal: Comput Biol Chem
DOI: 10.1016/j.compbiolchem.2020.107400
sha: fa00d0e8cfe90c1488676f7942b5a596647455f5
doc_id: 980321
cord_uid: lpjtdmij

A series of alkylated benzimidazole derivatives was synthesized and screened for their anti-HIV, anti-YFV, and broad-spectrum antiviral properties. The physicochemical parameters and drug-like properties of the compounds were assessed first, and then docking studies and MD simulations on HIV-RT allosteric sites were conducted to find the possible mode of their action. DFT analysis was also performed to confirm the nature of the hydrogen bonding interaction of active compounds. The in silico studies indicated that the molecules behaved like possible NNRTIs. The nature – polar or non-polar and position of the substituent present at fifth, sixth, and N-1 positions of the benzimidazole moiety played an important role in determining the antiviral properties of the compounds. Among the various compounds, 2-(5,6-dibromo-2-chloro-1H-benzimidazol-1-yl)ethan-1-ol (3a) showed anti-HIV activity with an appreciably low IC(50) value as 0.386 × 10(−5)µM. Similarly, compound 2b, 3-(2-chloro-5-nitro-1H-benzimidazol-1-yl) propan-1-ol, showed excellent inhibitory property against the yellow fever virus (YFV) with EC(50) value as 0.7824 × 10(−2)µM.

Most of the diseases caused by DNA and RNA viruses and the concomitant opportunistic infections represent a severe public health concern all over the world. Among viral diseases, AIDS (Acquired Immune Deficiency Syndrome) caused by human immunodeficiency virus (HIV) continues to be a deadly disease to date and stands as a major hazard to human health. Although several drugs have been discovered for the treatment of this dreadful disease, none of them were found to be efficient to kill the virus in toto, probably due to the emergence of drug resistance. To J o u r n a l P r e -p r o o f overcome this problem, several drug regimens are recommended during highly active antiretroviral therapy (HAART) used for the treatment of HIV/AIDS. As a result, improvement in adherence to medication is achieved and drug resistance has been lowered but, unfortunately, no regimen has been obtained to clean the virus out of patients completely (Tronchet and Seman, 2003; Martins et al., 2008; Zhan and Liu, 2011; Ouahrouch et al., 2016) . Thus, the emergence of drug resistance warrants the discovery of new drugs that standalone or as a drug regimen focused on a single target or different targets.

In the life cycle of HIV-1, there are several drug targets but reverse transcriptase (RT) is considered to be the most critical target for anti-HIV agents because it is the enzyme responsible for viral replication through the process of reverse transcription (Kumari and Singh, 2013; Tian et al., 2014; . HIV-RT is exposed to two major classes of inhibitors: nucleoside reverse transcriptase inhibitors (NRTIs) and non-nucleoside reverse transcriptase inhibitors (NNRTIs) (Meng et al., 2014; . NNRTIs bind to the flexible hydrophobic nonnucleoside inhibitor binding pocket (NNIBP) of RT allosteric site, 10Å away from the active site where NRTIs bind . Studies on mutation in RT have confirmed that some amino acid residues (Trp229, Phe227, Leu234 and Tyr318) present in the hydrophobic region of RT, play an important role in the accurate binding of NNRTIs to the protein (RT) and make the enzyme less vulnerable to mutation (Lansdon et al., 2010; La Regina et al., 2010) . So, a possible approach for designing new NNRTIs can be adopted to target this highly conserved region where the minimum possibility of mutation is expected.

Besides, the flexible and positional adaptable nature of NNRTIs plays an important role in enhancing the efficacy against a variety of drug-resistant strains of HIV-1 (Ferro et al., 2017) .

According to the drug design concept, a suitable inhibitor should have the efficiency to enter easily and occupy maximum space in the binding pocket to improve its potency.

The chemistry of benzimidazole moiety is of very high importance in the field of medicinal research (Digwal et al., 2016; El-Feky et al., 2015; Evans et al., 2015; Mavrova et al., 2013; Hranjec et al., 2012; El-Gohary and Shaaban, 2014; El-Gohary and Shaaban, 2015; Mahalakshmi and Chidambaranathan, 2015; Nara et al., 1989; Kus et al., 2008) . Several substituted benzimidazole derivatives (Monforte et al., 2010; Monforte et al., 2009; Monforte et al., 2008) have been reported as potent HIV-1 inhibitors (Fig. 1) where substitution at 5 th and 6 th position of J o u r n a l P r e -p r o o f benzimidazole moiety played an important role in antimicrobial/antiviral activity (Yadav and Ganguly, 2015; Akhtar et al., 2016; Francesconi et al., 2020; Cichero et al., 2017) .

We herein report the design of benzimidazole derivatives with a substituent at 5-, 6-position, and N-1 position, as possible NNRTIs against HIV-1. The synthesis of these molecules is reported in our previous work (Srivastava et al., 2018) as antibacterial agents and also given in supporting information. In the present article, we report the antiviral properties of benzimidazole derivatives.

The molecules were tested for their anti-HIV, anti-YFV and broad-spectrum antiviral properties against different DNA and RNA viruses, and the computational approaches, like molecular docking, DFT analysis, and molecular dynamics simulations were used to establish the mechanism of their mode of action.

All benzimidazole derivatives docked with HIV-RT protein interacted via hydrophilic interaction with Lys101, Lys103 or hydrophobic interaction with Tyr181, Tyr188, Phe227, Trp229 and Leu234 residues and this indicated positively towards their behaviour as NNRTIs. To compare the structural feature, docking, and biological results of these derivatives, nevirapine, a known NNRTI has been used as a reference drug.

General structures of nevirapine and N-1 alkylated benzimidazole derivatives are shown in Fig. 1.

The processes involved in the development of benzimidazole derivatives as NNRTIs against HIV-1 are summarized in Fig. 2. <Fig. 2>

Based on the Lipinski's rule of five (Lipinski et al., 2012) , twenty benzimidazole derivatives have been designed as NNRTIs against HIV-1 using Moleinspiration and Chemdraw software.

After performing the in silico structure-based experiments, eight (1(ab) -4(ab)) promising molecules have been selected for further studies. The physicochemical descriptors, like molecular weight (MW), hydrogen bond donor (HBD), hydrogen bond acceptor (HBA), the logarithmic partition coefficient in octanol/water, and logP (lipophilicity) have been analyzed to evaluate the drug-like properties of the molecules (Jacomini et al., 2016) . In silico ADMET (Adsorption, J o u r n a l P r e -p r o o f Distribution, Metabolism, Excretion, and Toxicity) predictions of all molecules are evaluated using DS 2.5 software to support the biological results.

Docking studies were performed using Discovery Studio 2.5 (DS 2.5, Accelrys Ltd., UK) and 3D X-ray crystal structure of RT protein (3MEC: PDB ID) was retrieved from Protein Data Bank (www.rcsb.org). Preparation, docking, and scoring of protein and molecular dynamics (MD) simulations of protein were carried out using the published procedure (Lansdon et al., 2010; Yadav et al., 2019) .

The target protein with docked ligand obtained from protein data bank was taken, extracted the ligand, and added the missing hydrogen atoms for preparing the receptor, using DS 2.5.

Optimization of the position of each atom was performed using the all-atom CHARMm forcefield with Adopted Basis set Newton Raphson (ABNR) minimization algorithm until the root mean square (r.m.s.) gradient for potential energy was <0.05 kcal/mol/Å (Brooks et al., 1983) . The target protein was minimized and defined as a receptor using the 'Binding Site' tool in DS 2.5. The space covered by the ligand in the receptor was defined as the binding site and an input site sphere having a radius of 5 Å was defined over the binding site. The minimized receptor created from the target protein was used for docking simulation.

The 3D structure of each ligand was built using the built-and-edit unit of DS 2.5 and then minimized using CHARMm forcefield by the ABNR method. A conformational search of the ligands was carried out by the MD approach. The ligand was heated up to 700K and then annealed up to 200K and this process was repeated thirty times. After thirty cycles, the conformation of each ligand was obtained and further subjected to local energy minimization using the ABNR method.

Each minimized conformation was then superimposed and the conformation with the lowest energy was selected as the most probable conformation.

To understand the behaviour and binding mode of benzimidazole derivatives, docking simulations were performed using DS 2.5 software. Docking is the computational prediction of the structures of ligand-protein complexes within a targeted binding site of protein receptors from the conformations of the unbound ligand and protein molecules (Kitchen et al., 2004) . The DS 2.5 J o u r n a l P r e -p r o o f Ligandfit protocol, combined with a shape comparison filter and a Monte Carlo transformational search, was used for docking of the ligand with targeted protein (Venkatachalam et al., 2003) .

Dreiding forcefield was used to refine the docked poses by rigid body minimization of the ligand for the grid-based calculated interaction energy (Mayo et al., 1990) . In the process of docking, the targeted protein receptor conformation was kept fixed. The minimization of the docked poses was carried out using the all-atom CHARMm forcefield and smart minimization method until the r.m.s. gradient for potential energy was <0.05 kcal/mol/Å. In the process of minimization, the atoms of ligand and the side chains of the residues of the receptor within 5 Å from the center of the binding site were kept flexible using simulation method (Molecular dynamics, energy minimization, and Monte Carlo simulation) (Yadav et al., 2019) .

Ligand scoring is an excellent approach to assess the true conformation and the binding tightness of ligand interacted into the active site of a target protein receptor (Kitchen et al., 2004) .

Scoring functions, like Interaction energy, Lig_Internal_Energy, Binding energy, Dock score, and Ludi_2 and Ludi_3 values were analyzed to find out the correct pose of protein-ligand complexes.

In docking simulation, Interaction energy is the energy of pairwise interaction of protein-ligand and complexes that explains the interaction of the ligand with the target protein. Ludi_2 and Ludi_3 scoring functions, binding energy, and EC50 values, respectively, of compounds were predicted (Kumar et al., 2010) . Docking simulations of all compounds into RT allosteric pocket generated ensembles of possible binding conformations. From these conformations, the most favourable conformation with the highest binding energy was selected for further analysis. The binding energy of ligands required for binding of ligands with the receptor in stable ligand-protein complexes (Lagos et al., 2008) was calculated by 'Calculate Binding Energy protocol' in DS 2.5 using the default settings (Bohm, 1994; Wang et al., 2003; Prathipati and Saxena, 2006) .

The validation of the docking protocol was performed by docking the native crystallized ligand (nevirapine) and molecules under study, in the allosteric binding site of RT (Reverse transcriptase) protein receptor. The comparison of docking results with the native crystallized ligand with root mean square deviation (RMSD) < 2 Å (Clark et al., 2016; Lam et al., 2018) , indicated that the used scoring function was appropriate. Docking values supported the hypothesis that experimental binding modes could be reproduced with correctness using this protocol.

MD simulations were carried out for the HIV-RT-compound 3a and HIV-RT-nevirapine complexes for a period of 50 ns using GROMACS (GROningen MAchine for Chemical Simulations) v5.1 molecular dynamics package (Pronk et al., 2013) .The unit cell defined as a cubical box, with a minimal distance of 10 Å from the protein surface to the edges of the box, was solvated using the Simple Point Charge (SPC) water model GROMOS96 53a6 force field (Oostenbrink et al., 2004) . Counter-ions were added to make every system electrically neutral at a salt concentration of 0.15 mol/L. Before the MD run, each system was subjected to energy minimization by employing the steepest descent integrator for 5000 steps with force convergence of <1000 kcal mol -1 nm -1 .

Thereafter, each protein-ligand complex was equilibrated for 5 ns using canonical (NVT) and isothermal-isobaric (NPT) ensembles. During equilibration, each system was coupled with the Berendsen temperature and Parrinello-Rahman pressure controllers, respectively, to maintain temperature 300 K and pressure 1 bar. The Particle Mesh Ewald (PME) algorithm (Essmann et al., 1995) was employed to deal with the long-range Coulomb interactions with a Fourier grid spacing of 0.12 nm. The short-range van der Waals interactions were assessed by the Lennard-Jones potential with a cut-off distance of 1 nm. All bond lengths were constrained by the linear constraint solver (LINCS) method (Hess et al., 1997) .

Subsequently, MD simulations were performed for 50 ns under micro-canonical ensemble by relaxing the couplings with the thermostats. In principle, the same protocol was applied to both systems. A time step of 2 fs was used and the coordinates were saved every 10 ps during the production run. For the structural analyses of every system, the resultant MD trajectories saved at the interval of 100 ps were analyzed using the built-in modules of GROMACS and visual molecular dynamics (VMD 1.9.1) (Schuler et al., 2001) . The 2-D plots depicting the intrinsic J o u r n a l P r e -p r o o f dynamical stabilities captured by the root mean square deviation (RMSD), and radius of gyration (Rg) of the complexes were generated by Grace 5.1.23 program.

MD simulations were performed on HIV-1 RT with compound 3a (showing anti-HIV activity) and reference drug nevirapine for investigating the interactions and quality of the binding and affinity with RT enzyme. Three independent MD simulations were used for this study. The MD results of nevirapine were used as a control.

Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) was taken into consideration for free energy calculation from MD trajectory (Aalten et al., 1995; Wang et al., 2018) . The binding free energy (ΔGbind) in a solvent medium was calculated as follows:

Where Gcomplex is the total free energy of the substrate-protein complex, Gprotein and Gligand are the total energies of protein and substrate alone in a solvent, respectively.

The free energies for each individual Gcomplex, Gprotein and Gligand were estimated by:

where p can be protein, ligand, or complex. EMM is the average molecular mechanics potential energy in vacuum and Gsolv is the solvation free energy.

The molecular mechanics potential energy was calculated in the vacuum as follows:

Where Ebonded or (Eint) is the total bonded interaction, which includes all bonded interactions like bond, angle, dihedral and improper interactions; Enon-bonded is the total non-bonded interaction consisting of both van der Waals (Evdw) and electrostatic (Eelec) interactions. Ebonded is always taken as zero.

The solvation free energy (Gsolv) was estimated as the sum of electrostatic solvation free energy (Gpolar) and nonpolar solvation free energy (Gnon-polar) as given below:

Where Gpolar, the polar solvation energy, was determined using the Poisson-Boltzmann (PB) linear equation and the non-polar contribution, Gnon-polar was estimated from the solvent-accessible surface area (SASA) as per the following equation:

Gnon-polar = γSASA + b

Where γ (coefficient related to surface tension of the solvent) = 0.02267 kJ/mol/Å 2 or 0.0054 kcal/mol/Å 2 and b = 3.849 kJ/mol or 0.916 kcal/mol

The binding free energies for all the complexes were calculated based on 5000 snapshots taken at an equal interval of time from 50 ns MD simulations.

The per-residue energy contribution was also computed to understand the contribution of individual amino acids to the total binding energy.

The molecular interactions of active compound 3a with HIV-RT protein through hydrogen bonding during docking and MD simulations were further confirmed by DFT optimized structurebased theoretical calculations (Sekhar et al., 2015; Parr and Yang, 1989) . The geometry optimization of complex RT:3a was carried out with the Gaussian 09 program (Frisch, et al., 2009 ).

The structure was optimized using the Density Functional Theory (DFT) method with the B3LYP/6-311g (d,p) basis set without any solvent to get the lowest energy structure. The energy minimized structures were further confirmed by the values of vibrational frequency. To confirm the stability of the complex RT:3a through hydrogen bonding NBO, QTAIM and NCI analyses were performed using the DFT method.

The natural bond orbital (NBO) analysis was carried out to investigate the stability of hydrogenbonded dimers. The oxygen/nitrogen atom having lone pair acts as a donor whereas X-H (X = N, O) acts as an acceptor in the strong intermolecular charge-transfer interaction. The second-order perturbation theory was used to calculate the stabilization energies of the dimers (Galembeck et al., 2014) .

Atoms in molecules (AIM), also cited as the quantum theory of atoms in molecules (QTAIM), analysis was performed using Multiwfn program (3.3.6 version) (Lu and Chen, 2012) to generate the bond path, a bond critical point (bcp) and its corresponding energies based on the quantum observable (electron density;  and the energies densities). The energies were calculated by the equation (Espinosa et al., 1998) as given below:

where, V(bcp) = Potential energy at bcp generated by the Multiwfn program.

The magnitudes of electron density () and the sign of Laplacian of electron density (∇ 2 ρ) for bcp of hydrogen bonds of interest in the complex of RT:3a were calculated using AIM calculations (Lakshmipriya et al., 2015) to confirm the presence of hydrogen bonding between compound 3a and Trp229 amino acid residue of HIV-RT protein.

The non-covalent interaction (NCI) analysis was used to detect the non-covalent interactions in real space based on the electron density and its derivatives and provided strong information about the steric repulsion, hydrogen bonding, and van der Waals interactions. It is a theoretical approach of the two-dimensional graphical plot of RDG (reduced density gradient) (s) versus sign (2) ×  (Johnson et al., 2010) .

Here, the Multiwfn program (Lu and Chen, 2012 ) was used to plot NCI, and VMD program (Humphrey et al., 1996) was used to plot colour filled isosurfaces for the complex RT:3a. In this analysis, the grid points were calculated and plotted for the two functions, i.e., sign (2) ×  (function 1) and RDG (s) (function 2). The function sign (2) ×  denoted the multiplication of electron density having the sign of second eigenvalue (2) of the electron density Hessian matrix.

The second function RDG was determined by the gradient of electron density (Johnson et al., 2010) . Percent inhibition was calculated by subtracting the above value from a hundred. Furthermore, the cytotoxicity assay was also performed to calculate the selectivity index of compounds.

To evaluate cytotoxic activity of test compounds on TZM-bl cells, an MTT [3-(4,5dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide; Sigma-Aldrich Inc.] assay was performed. The TZM-bl cells (3×10 4 cells per well) were seeded in 96-well plates (Greiner Bio-One, GmbH, Frickenhausen, Germany) in the presence of increasing concentrations of the test compounds and incubated in an incubator at 37°C with 5% CO2 for 48 h. the negative control included cells treated with solvent/medium. After incubation, cell viability was measured by the addition of 20 μL MTT (5 mg/mL per well) and incubated at 37°C for 3-5 h followed by the addition of MTT solvent (100 μL per well). The absorbance (OD) was read at 570 nm with the reference filter at 690 nm (Singh et al., 2016) . Experiments were performed three times, and the percent viability was calculated using the equation given below;

The compounds were also evaluated against different viruses such as herpes simplex virus-1 (HSV-1) (KOS strain), HSV-1 thymidine kinase-deficient, acyclovir-resistant (TK -Kos, ACV r ), 50% of the cell cultures) and the cell cultures were incubated in the presence of varying concentrations of the test compounds. Viral-induced cytopathicity was recorded as soon as it reached completion in the control virus-infected cell cultures that were not treated with the test compounds. Antiviral activity was expressed as the EC50 or compound concentration required to reduce the virus-induced cytopathicity by 50%.

Cytotoxicity of test compounds was expressed as the minimum cytotoxic concentration (MCC).

A colorimetric 3-(4,5-dimethythiazol-2-yl)-5-(3-carboxy-methoxyphenyl)-2-(4-sulfophenyl)-2Htetrazolium (MTS) assay system was used to calculate the 50% cytotoxic concentration (CC50) of test compounds responsible for 50% decrease in cell viability (Akula et al., 2017) .

All physicochemical parameters were found to be favourable for all benzimidazole derivatives (1(ab) -4(ab)). Data have been presented in Table 1 . From Table 1 , it was clear that all compounds were found within the acceptable range of Lipinski's rule of five and expected to behave as drug like molecules as nevirapine. All compounds had the potential to act as H-bond acceptors as well as donors due to the presence of -OH group, C-N and C-X dipole that helped to increase the solubility of compounds through ionization. In Table 1 , the lipophilicity of all compounds indicated that compounds had no problem to passage through the cell membrane. Lower TPSA values of all compounds proved that all compounds were capable to internalize well in cells. So, all compounds were expected to possess drug-like characteristics as shown in Table 1 .

To support the physicochemical properties of all benzimidazole derivatives (1(ab) -4(ab)), in silico ADMET parameters of all compounds were assessed and presented in Table 2 . Intestinal absorption was predicted to be more than 80% for all compounds like reference drug, which was greater than the acceptable range as shown in Table 2 . Compounds were considered either to behave as substrate or inhibitors of CYPsthe enzymes that play a vital role in the bioavailability of drugs. All compounds were found to be the inhibitors of CYP1A2, like the reference drug, as revealed in Table 2 .

All the ADMET parameters including absorption, distribution, metabolism, excretion, and toxicity were found to be permissible for all compounds and even better in some sense than that of reference drug (Table 2) .

J o u r n a l P r e -p r o o f <Table 2> (1(a -b) -4 (a -b) ) Substituted benzimidazoles 1-4 have been synthesized from 2-chloro-benzimidazole moiety.

The further modifications with different substituents at 5/6 position and N-1 position on substituted benzimidazoles are clearly outlined in Schemes 1 and 2. Purification and characterization data of final benzimidazole derivatives (1(ab) -4(ab) ) have been provided in supporting information.

All benzimidazole derivatives were docked into the NNIBP (non-nucleoside inhibitor binding pocket) of the X-ray crystal structure of RT (reverse transcriptase) protein to see their mode of interaction, like hydrogen bonding and non-covalent pi -pi, pi -cation interactions that stabilize the ligand-proteinin complexes (Chen et al., 2006; Gallivan and Dougherty, 1999) . The results obtained from docking simulations of all compounds with RT protein are shown in Table 3 .

From docking simulation, it was observed that the benzimidazole ring present in all compounds was placed inside the NNIBP of HIV-RT and surrounded by the hydrophobic residues, viz.

Val106, Val179, Tyr181, Tyr188, Leu234, phe227 and Trp229 and hydrophilic residues, viz.

Lys101 and Lys103. The hydrophilic substituent atoms/groups present on the benzimidazole core were situated towards the outer surface of receptor site marked with the hydrophilic amino acid residues. On the other hand the hydrophobic substituent atoms/groups present on the benzimidazole core were oriented towards the hydrophobic amino acid residues present inside the hydrophobic pocket of target RT protein. Binding mode and the hydrophobic or hydrophilic contacts of benzimidazole derivatives with particular amino acid residues of HIV-RT and the reference drug nevirapine into RT allosteric pocket are shown in Fig. 3 (a) and 3(b).

Scoring data of all molecules obtained from docking simulations are presented in Table 4 .

From (-42.48 Kcal/mol) showed interaction energy very close to that of the reference drug. Rest compounds showed lower interaction energy with respect to the standard drug used in molecular docking simulation.

The internal ligand energy (LIG) value of compounds 2b, 3a, 4a and 4b were somewhat comparable to that of the reference drug as shown in Table 4 .

Docking results illustrated that compounds 3a (-5.63 kcal/mol) and 4a (-5.60 kcal/mol) showed higher binding free energy than the reference (-5.57 kcal/mol) whereas compounds 3b (-5.37 Kcal/mol) and 4b (-5.45 Kcal/mol) showed binding energy very close to that of the reference drug, however, the other compounds showed a little lower values of binding energy than the reference drug.

From docking simulations on test compounds, it was established that higher value of binding energy was associated with a higher affinity towards stable protein-ligand complex, which, in turn, corresponded to lower EC50 value as presented in Table 4 .

Based on Dock Score, functional poses of ligand were evaluated and prioritized. It was confirmed from Table 4 that all compounds exhibited excellent docking scores ranging between 40.84 and 47.05 than the reference drug (40.60) except compounds 1a and 1b. All molecules had a similar pattern of binding into the RT allosteric site and also had a significant binding affinity with RT core, just like the reference drug. Therefore, the molecules were studied further to explore their structure-activity relationships and inhibitory activity as NNRTIs against HIV-1.

All N-1 alkylated benzimidazole derivatives were screened for their ability to inhibit the HIV-1 and the results obtained are summarized in Table 5 .

From Table 5 , it was observed that amongst compounds 1(ab) -4(ab) only compound 3a with bromo groups at 5 th and 6 th position and -CH2CH2OH group at N-1 position showed significant inhibition ( IC50 = 0.386×10 -5 µM) against HIV-1 as compared to other compounds having -H, -NO2 substituents at 5 th , 6 th position and -CH2CH2CH2OH at N-1 position. However, the selectivity index (SI = 2.73) of this compound was much lower than the reference drug. All compounds showed high to moderate cytotoxicity against TZM-bl cells (Table 5 ). Fig. 4 showed the percentage inhibition of HIV-1 infection and cytotoxicity of compound 3a at the tested concentration.

Di-bromo groups rendered higher hydrophobicity than nitro or di-nitro groups; therefore, compound 3a had higher chances for hydrophobic interactions and higher chances for H-bonding due to the presence of polar group at N-1 position with HIV RT and this could be the reason for its anti-HIV activity. However, the absence of pi-pipi interaction with Tyr181 and Tyr188 amino acid residues of the hydrophobic pocket of RT protein might be the reason for lowering its activity.

SAR studies (Fig. 5) confirmed the significance of the substituents at fifth, sixth and N-1 position of benzimidazole nucleus for anti-HIV activity. The structure-activity relationship study revealed that halogen groups at positions 5 th , 6 th and -CH2CH2OH group at position N-1 in 2-Cl-1H-benzimidazole nucleus favoured the inhibitory activity as compared to other substituents at 5 th , 6 th position (-H, -NO2 group) and N-1 position (-CH2CH2CH2OH group).

The broad-spectrum antiviral activities of N-1 alkylated benzimidazole derivatives were tested against different DNA/RNA viruses and the results are given in Table 6 .

From Table 6 , it was confirmed that all test compounds failed to show any antiviral activity against different viruses except compound 2b, which showed an inhibitory property against YFV in Vero cell cultures (Table 6 ). However, no activity was found in HEL, HeLa and MDCK cells.

It can be observed that compound 2b having -NO2 group at 5 th position and -CH2CH2CH2OH group at N-1 position of benzimidazole nucleus had an ability to inhibit the replication of YFV at micromolar concentrations. Compound 2b showed YFV inhibition at EC50 values ranging from 0.0234µM to 0.7824×10 -2 µM and SI values ranging from 16.64 to 50, which was higher than the reference drug ribavirin but lower than the DS-10.000 and mycophenolic acid.

To investigate the stability of RT:nevirapine and RT:3a complexes, the RMSD values of these complexes were calculated (Fig. 6) . Apart from assessing the equilibration, quality of the run and convergence of MD trajectories, RMSD is also useful in investigating the stability of a protein in complex.

From Fig. 6 , it was clear that there was a fluctuation in the RMSD value of RT:3a complex initially up to 1.6 nm until 15 ns, and then it got converged and remained stable throughout the 50 J o u r n a l P r e -p r o o f ns duration of MD simulation, which indicated the attainment of the stability of RT:3a complex after 15 ns.

In case of RT:nevirapine complex, fluctuation of RMSD could be seen at 20 ns, 25 ns and 35 ns, after that it remained stable till 50 ns.

The radius of gyration (Rg) describes the level of compaction of protein. It is defined as the mass-weighted root-mean-square distance for a collection of atoms from their common center of mass. Hence, the trajectory analysis for Rg provided conformational variation, compactness and tertiary structural volume of the protein-ligand complexes. During 50 ns MD simulations, the radius of gyration for RT:3a and RT:nevirapine complexes was calculated and analyzed to see the effects of the binding mode of ligand on structural features of RT backbone (Fig. 7) .

Like RMSD, initial fluctuation of Rg till 15 ns could be seen in the case of RT:3a complex, after that it became constant and remained stable throughout the 50 ns MD simulation. Similarly, in case of RT:nevirapine complex, fluctuation could be seen during 18 to 35 ns of MD simulation time period and remained stable after that.

Trajectory analysis, after 50 ns, indicated that during simulations both compounds maintained their interactions with RT as observed in their docked poses (Fig.8) .

A combination of MM/PBSA and MD method was used to assess the binding free energy of RT:nevirapine and RT:3a complexes to corroborate the results. The snapshots were extracted from the 50 ns of MD simulation of each component to analyze the binding free energy of both complexes. All the energy terms of MD simulations are listed in Table 7 & Fig. 9 .

In Table 7 , the binding free energy (∆Gbinding) of RT:nevirapine and RT:3a complexes was -85.31 and -88.33 kJ/mol, respectively. From Table 7 , it was observed that compound 3a showed almost similar calculated binding energy as that of the reference drug nevirapine. The calculated van der Waals contributions (∆Evdw) to the binding free energy in the RT:3a (-173.32 kJ/mol) complex was more than that of RT:nevirapine complex (-162.89 kJ/mol). A sum of the calculated solvation energy (∆Esolv), van der Waals energy (∆Evdw) and electrostatic energy (∆Eelect) favoured J o u r n a l P r e -p r o o f the binding of compound 3a to RT allosteric site and also indicated the possibility of this molecule working as NNRTI.

In addition to the information presented in Fig. 9 & Table 7 , the contributions of individual amino acids to the binding free energy (∆Gbind) were also computed using the MM/PBSA method and presented in Fig. 10 . For clarity, only actively contributing residues towards the positive and negative binding energies are presented.

The per-residue interactional energy profiles revealed that Leu100, Lys101, Val106, Tyr181, Tyr188, Val189, Phe227, Trp229, Leu234 and Tyr318, actively participated in interaction and resulted in a stronger binding and stability in both the complexes. Interestingly, these residues were also interacting during the docking analysis, validating the molecular docking result.

The contribution of more than 5 kJ/mol for both RT:3a and RT:nevirapine complexes by two active residues, viz. Leu100 and Tyr181, indicated the importance of these residues in imparting stability to the complexes formed.

The stability of complex (RT:3a), through hydrogen bonding between compound 3a and Trp229 amino acid residue of HIV-RT protein in docking simulation, was further confirmed by NBO analysis. In NBO analysis of RT:3a complex, an important interaction from n*, i.e., (N(13) H(43)-O(34)) was observed. The n* interaction of a hydrogen bond involves the intermolecular delocalization (charge transfer) between a lone pair donor n and an antibonding acceptor * of an adjacent molecule, which stabilizes the structure of the protein-ligand complex.

From NBO analysis, it was found that the presence of n* interaction of hydrogen bonding was due to the transfer of lone pair electron from N-atom of Trp229 of HIV-RT to the antibonding orbital (*) of O-H bond that stabilized the protein-ligand complex (RT:3a). The stabilization energy due to the process of electron transfer for the hydrogen bond is given in Table 8 complex. The bcp and bond path of intermolecular hydrogen bond for the RT:3a complex is represented in Fig. 12 . The calculated values of electron density () and Laplacian of electron density (∇ 2 ρ) at each intermolecular bcp and the calculated value of hydrogen bond interaction energy (EHB) for the located bcp of hydrogen bond interest are expressed in Table 9 .

The AIMS analysis revealed the electron density value well in the expected range (0.0102-0.0642 a.u.) for a hydrogen bond (Shahi and Arunan, 2014) . The positive value of Laplacian of electron density at bcp showed that the type of interaction between compound 3a and Trp229 amino acid residue of HIV-RT was hydrogen bonding, which stabilized the formation of a proteinligand complex.

The NCI analysis was performed to affirm the presence of the hydrogen bond in the complex RT:3a. The functions carrying information about the presence of hydrogen bonding, the plot of RDG (s) versus sign (2) × , and the coloured gradient isosurfaces for the complex RT:3a are represented in Fig. 13(a) and 13(b) , respectively.

The one spike or negative peak with sign (2) ×  ≈ -0.020 on the left-hand side in Fig. 10(a) represented one hydrogen bond, i.e., (N H-O) in the complex RT:3a and it can be seen in Fig.   10 (b) as isosurface of blue colour. In Fig. 10(b) , the red and green colour isosurfaces showed the steric interactions and van der walls interactions, respectively. N-1 alkylated benzimidazole derivatives 1(ab) -4(ab) bearing H/NO2/Br groups at fifth and sixth position and hydroxyl alkyl groups at N-1 position was synthesized based on in silico studies. Molecular docking results indicated that binding mode of these compounds to HIV-RT allosteric site was similar to that of the reference drug, and hence, these compounds were supposed to behave as NNRTIs on the pattern of nevirapine. The compounds were screened against HIV-1 in TZM-bl cells but the expected inhibitory activity was not observed. Only the compound 3a showed anti-HIV activity with the IC50 value of 0.386×10 -5 µM but the lower SI value of this J o u r n a l P r e -p r o o f compound marred its importance. MD simulations results and DFT analysis indicated that compound 3a formed a stable complex with HIV RT, like the reference drug nevirapine. These results provided new insights to modify benzimidazole derivatives for further development of potential anti-HIV agents.

The broad-spectrum antiviral screening in Vero cell cultures revealed that compound 2b had excellent inhibition properties against the yellow fever virus with EC50 value 0.7824×10 -2 µM, equivalent to an SI value of 50.

The results found in the present study paved the way for further improvement in the antiviral property of the molecules through computational simulations.

Ritika Srivastava synthesized and purified all the compounds along with drafted the manuscript. 

The authors have no conflict of interest.

General structure of nevirapine and new designed benzimidazole derivatives.

Flow chart of processes involved in the development of benzimidazole derivatives.

Docked poses with key contacts and expanded views of molecular docking interactions of compounds 1a, 1b, 2a, 2b, 3a, 4a and 4b within NNIBP of RT protein.

Surface view of compound 3a and reference drug nevirapine within NNIBP of RT protein.

Anti-HIV activity and cytotoxicity of compound 3a at different concentrations used in assay.

SAR derived from screening of test compounds against HIV-1.

RMSD plots of the HIV-RT complexes of nevirapine and compound 3a, (shown in black and red, respectively), as a function of time.

Rg (radius of gyration) plots of RT:nevirapine (black) and RT:compound 3a (red) complexes, as a function of time. 1a, 1b, 2a, 2b, 3a, 4a , and 4b within NNIBP of RT protein. 

Reference drug (nevirapine) Fig. 8(a) . The binding mode of compound 3a and reference drug in hydrophobic pocket of RT allosteric site after 50 ns MD simulation. (Hydrogen bond and pi-pi interactions are shown by the green colour line in compound 3a and reference drug, respectively)

Compound 3a with nevirapine Fig. 8(b) . Superimposition of compound 3a at different time intervals; 10, 20 and 50 ns during MD simulations and its superimposed image with nevirapine after 50 ns.

RT:nevira pine J o u r n a l P r e -p r o o f Table 1 Physiochemical data of compounds 1(ab) -4(ab).

In Silico ADMET properties of compounds 1(ab) -4(ab). Table 3 Docking interaction of compounds 1(a -b) -4(a -b) within NNIBP of RT protein in ligand -receptor complex.

Docking scores of compounds 1(a -b) -4(a -b) within NNIBP of RT protein in ligand -receptor complex.

Anti-HIV-1 activity of compounds 1(a -b) -4(a -b) evaluated in TZM-b1 Cell cultures.

Antiviral activity of compounds 1(a -b) -4(a -b) evaluated in Vero cell cultures.

Molecular energy data of RT:nevirapine and RT:3a complexes.

Natural bond orbital analysis of orbital overlap hydrogen bonding interaction between the Trp 229 residue of RT and compound 3a.

Parameters of the bcp observed between the Trp 229 residue of RT and compound 3a: distance, hydrogen bond interaction energy, electron density and the Laplacian of electron density. Table 1 Physiochemical data of compounds 1(a -b) -4(a -b). Table 2 In Silico ADMET properties of compounds 1(a -b) -4(a -b).

Metabolism Excretion Toxicity Water solubility (log mol/L) Table 3 Docking interaction of compounds 1(a -b) -4(a -b) within NNIBP of RT protein in ligand -receptor complex. 

Lys 103 Trp 229

Lys 103: N -2b: O12 Trp 229: N -2b: O17 

Tyr 181 (' π-π') Tyr 188 (' π-π')

Tyr 181 -ref.

Tyr 188 -ref. Table 4 Docking scores of compounds 1(a -b) -4(a -b) within NNIBP of RT protein in ligand -receptor complex. 0.27×10 -5 I.E. = Interaction energy, DS = Dock score, Ludi_2 and Ludi_3 = empirical scoring functions, ∆G = Binding energy, LIG = Lig-Internal_Energy, Predicted EC50 = Predicted 50% effective concentration of given compounds essential to reduce HIV-1 replication (μM)

Anti-HIV-1 activity of compounds 1(a -b) -4(a -b) evaluated in TZM-b1 Cell culture. IC50 = 50% inhibitory concentration, which is the concentration required in 50% inhibition in HIV infection. CC50 = 50% cytotoxic concentration, which is the concentration required to reduce TZM-bl cell viability by 50%. SI = Selectivity index ratio CC50/IC50. Table 8 Natural bond orbital analysis of orbital overlap between the Trp 229 residue of RT and compound 3a suggesting hydrogen bonding.

Parameters of the bcp observed between the Trp 229 residue of RT and compound 3a: distance, hydrogen bond interaction energy, electron density and the Laplacian of electron density. 

Therapeutic evolution of benzimidazole derivatives in the last quinquennial period

Facile functionalization at the C4 position of pyrimidine nucleosides via amide group activation with (benzitriazol-1-yloxy)tris(dimethyamino)phosphonium hexafluorophosphate (BOP) and biological evaluation of the products

The development of a simple empirical scoring function to estimate the binding constant for a protein-ligand complex of known three-dimensional structure

CHARMM: a program for macromolecular energy, minimization, and dynamics calculations

On evaluating molecular-docking methods for pose prediction and enrichment factors

Acylthiocarbamates as non nucleoside HIV-1 reverse transcriptase inhibitors: docking studies and ligand-based CoMFA and CoMSIA analyses

Computational studies of the binding mode and 3D-QSAR analyses of symmetric formimidoester disulfides: a new class of non-nucleoside HIV-1 reverse transcriptase inhibitor

Thiocarbamates as non nucleoside HIV-1 reverse transcriptase inhibitors: Docking-based CoMFA and CoMSIA analyses

Benzimidazole-based derivatives as privileged scaffold developed for the treatment of the RSV infection: a computational study exploring the potency and cytotoxicity profiles

Prediction of protein-ligand binding poses via a combination of induced fit docking and metadynamics simulations

VOSO4 catalyzed highly efficient synthesis of benzimidazoles, benzothiazoles, and quinoxalines

Synthesis, molecular docking and anti-inflammatory screening of novel quinoline incorporated pyrazole derivatives using the Pfitzinger reaction II

Antimicrobial and antiquorum-sensing studies. Part 2: synthesis, antimicrobial, antiquorum-sensing and cytotoxic activities of new series of fused

Antimicrobial and antiquorum-sensing studies. Part 3: synthesis and biological evaluation of new series of

Hydrogen bond strengths revealed by topological analyses of experimentally observed electron densities

A smooth particle mesh Ewald method

Benzimidazole analogs inhibit respiratory syncytial virus G protein function

Searching for novel N1-substituted benzimidazol-2-ones as non-nucleoside HIV-1 RT inhibitors

Synthesis and Biological Evaluation of Novel (thio)semicarbazone-Based Benzimidazoles as Antiviral Agents against Human Respiratory Viruses

Effects of the protonation state in the interaction of an HIV-1 reverse transcriptase (RT) amino acid, Lys101, and a non nucleoside RT inhibitor, GW420867X

Cation-pi interactions in structural biology

LINCS: A linear constraint solver for molecular simulations

Synthesis, crystal structure determination and antiproliferative activity of novel 2-amino-4-aryl-4,10-dihydro

VMD: Visual molecular dynamics

Synthesis and evaluation against Leishmania amazonensis of novel pyrazolo[3,4-d]pyridazinone-N-acylhydrazone-(bi)thiophene hybrids

Revealing Non-Covalent Interactions

Docking and scoring in virtual screening for drug discovery: Methods and applications

In silico structure-based design of a novel class of potent and selective small peptide inhibitor of Mycobacterium tuberculosis dihydrofolate reductase, a potential target for anti-TB drug discovery

Anti-HIV drug development: Structural features and limitations of present day drugs and future challenges in the Successful HIV/AIDS treatment

Synthesis and antioxidant properties of novel N-methyl-134-thiadiazol-2-amine and 4-methyl-2H-124-triazole-3(4H)-thione derivatives of benzimidazole class

Looking for an active conformation of the future HIV type-1 non-nucleoside reverse transcriptase inhibitors

Docking and quantitative structure-activity relationship studies for the bisphenylbenzimidazole family of non-nucleoside inhibitors of HIV-1 reverse transcriptase

Three centered hydrogen bonds of the type C=O…H(N)…X-C in diphenyloxamide derivatives involving halogens and a rotating CF3 group: NMR, QTAIM, NCI and NBO studies

Ligand-biased ensemble receptor docking (LigBEnD): a hybrid ligand/receptor structure-based approach

Crystal structureof HIV-1 reverse transcriptase with etravirine (TMC 125) and relpivirine (TMC 278): Implications for drug design

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

Multiwfn: A multifunctional wave function analyzer

Synthesis, spectral characterization and antimicrobial studies of novel benzimidazole derivatives

The current status of the NNRTI family of antiretrovirals used against HIV infection

Design, synthesis and antiproliferative properties of some new 5-substituted-2-iminobenzimidazole derivatives

DREIDING: a generic force field for molecular simulations

Design and synthesis of a new series of modified CH-diarylpyrimidines as drugresistant HIV non-nucleoside reverse transcriptase inhibitors

Novel N1-substituted 1,3-dihydro-2H-benzimidazol-2-ones as potent non-nucleoside reverse transcriptase inhibitors

Novel 1,3-dihydro-benzimidazol-2-ones and their analogues as potent non-nucleoside HIV-1 reverse transcriptase inhibitors

Design, synthesis, and structure-activity relationships of 1,3-dihydrobenzimidazol-2-one analogues as anti-HIV agents

-Pyridylaminomethyl)benzimidazole derivatives having antiviral activity. US patent

A biomolecular force field based on the free enthalpy of hydration and solvation: The GROMOS force-field parameter sets 53A5 and 53A6

Design, synthesis and antiviral activity of novel ribonucleosides of 1,2,3-Triazolylbenzyl-aminophosphonates

Density-Functional Theory of Atoms and Molecules

Evaluation of binary QSAR models derived from LUDI and MOE scoring functions for structure based virtual screening

GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation tool kit

An improved GROMOS96 force field for aliphatic hydrocarbons in the condensed phase

Density functional theory, natural Bond Orbital and atoms in molecule analyses on the hydrogen bonding interactions in 2-chloroaniline-carboxylic acids

Hydrogen bonding, halogen bonding and lithium bonding: an atoms in molecules and natural bond orbital perspective towards conservation of total bond order, inter-and intra-molecular bonding

Design and anti-HIV activity of aryl sulphonamides as non-nucleoside reverse transcriptase inhibitors

Synthesis, antibacterial activity, synergistic effect, cytotoxicity, docking and molecular dynamics of benzimidazole analogues

Used heterocyclic compounds bearing bridgehead nitrogen as potent HIV-1 NNRTIs. Part 1: Design, synthesis and biological evaluation of novel 5,7-disubstituted pyrazolo[1,5-a]pyrimidine derivatives

Non nucleoside inhibitors of HIV-1 reverse transcriptase: from the biology of reverse transcriptase to molecular design

Essential dynamics of the cellular retinol-binding protein evidence for ligand-induced conformational changes

LigandFit: a novel method for the shape-directed rapid docking of ligands to protein active sites

Recent Developments and Applications of the MMPBSA Method

Comparative evaluation of 11 scoring functions for molecular docking

Structure activity relationship (SAR) study of benzimidazole scaffold for different biological activities; A mini-review

silico studies on new oxathiadiazoles as potential anti-HIV agents

Novel HIV-1 non-nucleoside reverse transcriptase inhibitors: a patent review

The author Ritika Srivastava acknowledges financial support in the form of Research 

Yes Water Solubility = < -4 soluble; Intestinal absorption = Below 30% indicates poor absorbance; Blood brain barrier Permeability = < -1considered poorly distributed to the brain; CNS (Central Nervous System) permeability = > -2 considered to penetrate the CNS; Total Clearance (logCLtot) = Lower value indicates high drug half lifetime; LOAEL (Lowest Observed Adverse Effect) = Lower value predicts minimum toxicity.