Abstract
Hepatitis C virus (HCV) is responsible for various clinical conditions ranging from acute viral hepatitis to chronic liver disease and cirrhosis, causing liver cancer. The inhibition of HCV NS5B polymerase is a major step towards design of anti-HCV drugs. This study is devoted to two potential inhibitors of NS5B polymerase using quantum chemical and molecular docking methods. The structures of these molecules have been optimized using density functional theory at B3LYP/6–311++G (d,p) level. The vibrational spectral analyses of these two molecules have been performed in detail. The calculated vibrational spectral studies of both the title compounds are in good agreement with previous result. The frontier orbital surfaces and molecular electrostatic potential surfaces have been analyzed. The energy gap between HOMO and LUMO suggest that the compd2 is more reactive than compd1. The molecular electrostatic potential plots suggest the suitable nucleophilic and electrophilic sites in title compounds. Various electronic and thermodynamic parameters of both molecules have been calculated and compared. Finally, the inhibition activity of these compounds with the help of molecular docking studies, against NS5B polymerase enzyme has been explored. The result of docking studies suggests that these compounds are emerging as an effective potent inhibitor of HCV NS5B polymerase and can be useful in the design of novel anti-HCV drugs.
Introduction
Hepatitis C virus (HCV) has been a serious medical issue, which results in several clinical conditions varying from acute viral hepatitis to chronic liver disease and cirrhosis. It is the major factor responsible for liver cancer and around 3% of the world’s population, are chronically affected with HCV infection such that 3–4 million people are newly infected each year [1]. Although there is no vaccine available to treat HCV infection, the current clinical method of care includes a combination of pegylated interferon (peg-IFN) with the nucleoside analog ribavirin (RBV) [2]. It should, however, be noted that the identification of more effective HCV inhibitors is still required due to low response rates of peg-IFN/RBV therapies and their significant side-effects, particularly for patients infected with genotype 1 HCV [3]. It has been reported [4] that the activity of the virally encoded HCV NS5B polymerase is essential for HCV replication making it an attractive target for the development of novel anti-HCV drugs. In particular, two class of compounds namely, monocyclic 5,6-dihydro-1H-pyridin-2-ones [5] and 1,1-dioxo-1,4-dihydro-1λ6-benzo[1,2,4]thiadiazine moiety linked to a bicyclic 5,6-dihydro-1H-pyridin-2-one [6] have been reported to be potent inhibitors of NS5B polymerase. In this study, we will explore the chemical properties, reactivity and inhibition activity of two compounds belonging to latter class (as described above) against NS5B polymerase enzyme using quantum chemical and molecular docking calculations.
Quantum chemical methods have already proven to be very useful in determining the molecular structure as well as elucidating the electronic structure and reactivity of a variety of biologically active molecules [7–10]. Although, there are several quantum chemical methods available nowadays, densities functional theory (DFT) based methods have gained much popularity due to their better compromise between computational cost and accuracy. DFT methods generally provide the electronic energy and related parameters as a function of electron density, which is itself a function of position coordinates. These methods differ in the way they approximate the exchange and correlation parts of the energy functional. Currently available DFT methods are fully capable to reproduce experimental data such as structural parameters (bond lengths, bond angles, etc.) and spectral parameters (FT-IR, NMR, UV-vis, etc.) and to provide reasonable predictions of the same in the absence of experimental data. Considering these facts, we present DFT calculations on chemical properties of two NS5B polymerase inhibitors belonging to 1,1-dioxo-1,4-dihydro-1λ6-benzo[1,2,4]thiadiazine linked with a bicyclic 5,6-dihydro-1H-pyridin-2-one series of compounds.
Computational methods
In the present study, we have performed all the density functional theory based calculations using the Becke’s 3-parameter hybrid functional for exchange part and the Lee-Yang-Parr functional for correlation (B3LYP) [11, 12] with 6–311++G(d,p) basis set using the Gaussian 09 suite of program [13]. The basis set 6–311++G(d,p) is a triple-split valance basis set that increases the flexibility of the valance electrons. It is useful for most of the studies involving medium-size system [14]. The vibrational frequencies are also calculated at the same level of theory and tabulated in Table 2(a, b). The vibrational wavenumber assignments have been carried out by combining the result of the Gauss view 05 program [15] symmetry considerations and the VEDA04 program [16].
Results and discussion
Optimized geometry
The optimized geometries of title compounds are calculated at B3LYP/6–311++G(d,p) level of theory are shown in Fig. 1(a) and 1(b) with labeled atoms. Two compounds differ only in left moieties attached to six-membered ring containing one N and two substituted O atoms. In compound 1 (compd1), there is one fluorinated methylbenzene attached to N and ethyl moiety attached to neighboring C atom as shown in Fig. 1(a). In compound 2 (compd2), on the contrary, five-membered ring is attached to N and alkyl chain is attached to C atom, as evident from Fig. 1(b). The selected bond-lengths and bond-angles calculated for equilibrium geometries of compd1 and compd2 are listed in Table 1. From Table 1, it is clear that most of the bond lengths and bond angles are similar in compd1 and compd2.
(a) Optimized geometry of Compd1 with labeled atom. (b) Optimized geometry of Compd2 with labeled atom.
Selected bond lengths (Angstroms), bond angles (degrees) of compd1 and compd2, calculated at the B3LYP/6–311++G(d,p) level
The vibrational assignments for compd1 and compd2 along with potential energy distribution (PED%) are presented in Tables 2 and 3, respectively. The Calculated IR spectra of title compounds are shown in Fig. 2. Below, we discuss and analyse the significant modes of vibration in the molecules under study.
Vibrational analysis of some selected modes of Compd1 at B3LYP/6–311++G(d,p) level
Vibrational analysis of some selected modes of Compd1 at B3LYP/6–311++G(d,p) level
Vibrational analysis of some selected modes of Compd2 at B3LYP/6–311++G(d,p) level
*PED < 15% have been neglected in assignments. Abbreviations used: R- ring, ν- stretching, νs- symmetric stretching, νas- antisymmetric stretching, σ- scissoring, ρ- rocking, τ- twisting, τi- in-plane torsion, τo- out-of-plane torsion.
Calculated vibrational IR spectra of Compd1 and Compd2 at B3LYP/6–311++G(d,p) level.
The N–H stretching vibrations are normally viewed in the region 3300–3600 cm–1 [17]. For compd1, the N-H stretching vibration observed at 3455 cm–1 and 3163 cm–1 while 3458 cm–1 and 3230 cm–1 is observed for compd2. The other modes of vibration such as rocking, scissoring for compd1 and compd2 are observed at lower frequency below 1500 cm–1.
C–H bands usually present in the region 2800–3200 cm–1 [18, 19]. In the present study, strong C–H stretching vibrations for both the compounds are observed in the range 2880–3120 cm–1. All theoretical wavenumbers due to C–H stretching by B3LYP/6–311++G(d,p) are in good agreement with the results of Sawant and coworkers [20]. Mixing of C–H and C–N rocking and scissoring vibrations are observed in the middle region of the spectra.
C = O and S = O vibrations
In literature, strong C = O bands are found in the region 1800–1690 cm–1 [21, 22]. The calculated frequency of the C = O stretching reduces due to π–π bond formed between C and O. The nature of the carbonyl group is determined by the lone pair of electrons on oxygen. In the present study, for compd1, a very strong C = O band is shown at 1672 cm–1 and for compd2 at 1569 cm–1. The S = O stretching vibrations are observed below 1500 cm–1 for both the compounds.
Ring vibrations
The C–C stretching vibrations give rise to characteristics bands in the spectral range 1600–1400 cm–1. Normally the C–C stretching vibrations are expected within the middle region [23, 24]. For the compd1 and compd2, the C–C stretching vibrations in aromatic compounds are in the range 1472–1601 cm–1. Most of the ring vibrational modes are affected by the substitution in the aromatic ring. The in-plane and out-of-plane modes of vibration are seen in the lower region below 1000 cm–1.
HOMO-LUMO and MESP surfaces
The highest occupied molecular orbital (HOMO), lowest unoccupied molecular orbital (LUMO) and their properties such as energy are very useful for physicist and chemists and are very important parameters for quantum chemistry. This is also useful for predicting the most reactive position in π-electron system and also explains several types of reaction in conjugated system [25]. The conjugated molecules are characterized by a small highest occupied molecular orbital and lowest unoccupied molecular orbital (HOMO-LUMO) separation, which is the result of significant degree of intra-molecular charge transfer from the end-capping electron-donor groups to the efficient electron acceptor groups through π-conjugated path [26]. The HOMO-LUMO energy gaps are used to prove the bioactivity from intra-molecule charge transfer [27, 28]. The HOMO and LUMO plots of compd1 and compd2 have also been plotted with the help of GaussView 5.0 program and are shown in Fig. 3. The energy gap of compd1 is 4.41 eV and that of compd2 is 4.34 eV. This may suggest that compd2 is slightly more chemically reactive than compd1.
HOMO-LUMO plots of (a) Compd1 (iso value = 0.02 a.u.) (b) Compd2. (iso value = 0.02 a.u.).
The importance of molecular electrostatic potential (MESP) lies in the fact that it simultaneously displays molecular size, shape as well as positive, negative and neutral electrostatic potential regions in terms of color grading, which is very useful in the investigation of the most probable binding receptor site along with the size and shape of the molecule [29, 30]. The MESP plots of compd1 and compd2 are displayed in Fig. 4. In color grading scheme of MESP, the blue represent the most electropositive, i.e., electron poor region, whereas the red corresponds to the most electronegative centre, i.e., electron rich region. From the Fig. 4, it is clear that the most electronegative region for both the compounds is located over the oxygen atom doubly bonded with sulphur atom, indicating a possible site for electrophillic attack.
MESP plots of title compounds, (a) Compd 1, (b) Compd 2.
We have calculated various electronic parameters of title compounds, viz. ionization potentials (I), electron affinity (A), absolute electro negativity (χ) and chemical hardness (η) etc. at B3LYP/6–311++G(d,p) level. I and A are calculated as the negative of energy eigen-values of HOMO and LUMO respectively. χ and η can be calculated by using finite–difference approximations [31] as η= ½ (I–A) and χ= ½ (I + A). These parameters often used to describe chemical reactivity of molecules are listed in Table 4. It is clear from the Table 4 that compd1 possesses slightly higher values as compared to compd2 and hence high electronegativity. The thermodynamic parameters are also listed in Table 4, i.e., zero point energy (ZPE), thermal energy at room temperature (E), heat capacity (Cv) and entropy (S) for both molecules at the same level of theory. These parameters are related to one another via standard thermodynamical relations and can be very useful in estimating reaction paths of molecules.
Electronic and thermodynamic parameters calculated at B3LYP/6–311++G (d,p) level
Electronic and thermodynamic parameters calculated at B3LYP/6–311++G (d,p) level
It has been reported that molecules containing a 1,1-dioxo-1,4-dihydro-1-benzo[1,2,4] thiadiazine moiety linked to a bicyclic 5,6-dihydro-1H-pyridin-2-one are potent inhibitors of the NS5B polymerase [32]. We have performed the molecular docking simulation of the compounds with HCV NS5B inhibitor receptor by Hex program (version 8.0.0) [33]. The 3D crystal structure of target protein was obtained from Protein Data Bank (PDB ID: 3HHK) [34]. Hex can calculate protein-ligand docking, assuming the legend is rigid, and it can superimpose pairs of molecule using only knowledge of their 3D shapes. The docking score can be approximated to an interaction energy value (e-value), which we seek to minimize. The more negative the e-value, the more efficient will be the docking process. Note that it is merely a model which may provide the binding affinity of a particular site in terms of e-value. The docking of title compounds are displayed in Figs. 5 and 6. The total e-value obtained is –330.34 a.u. for compd1 and –374.01 a.u. for the compd2. From this result we may conclude that compd2 is more effective than compd2 as an inhibitor of HCV NS5B polymerase.
Docking of compd1 into 3HHK inhibitor. The ligand compd1 is displayed by ball and stick model in carton view (a) and by solid surface (b).
Docking of compd2 into 3HHK inhibitor. The ligand compd2 is displayed by ball and stick model in carton view (a) and by solid surface (b).
This paper presents a comprehensive computational study on two HCV NS5B polymerase inhibitors. All calculated wavenumbers are real in nature for both the molecules thus both compounds are stable. The detailed descriptions of the vibrational spectra of these molecules have been done with the help of normal modes analysis. Reactivity plays a key role in specific chemical reactions. The lower value of frontier orbital energy gap in case of compd1 suggests a more reactive nature as compared to compd2. We have also discussed sites for both molecules during electrophilic, nucleophilic, and radical attacks with the help of global reactivity descriptors. The result of molecular docking studies suggests that these compounds are emerging as an effective potent inhibitor of HCV NS5B polymerase. The findings from this study may provide a suitable path for researchers in future.
Footnotes
Acknowledgments
GT is thankful to Dr. R. P. Ojha, Department of Physics, D.D.U. Gorakhpur University for helpful discussion and needful guidance. AKS is thankful to Science and Engineering Research Board (SERB) for National Postdoctoral Fellowship [Grant No. PDF/2016/001784].