QSAR modeling, molecular docking, dynamic simulation and ADMET study of novel tetrahydronaphthalene derivatives as potent antitubercular agents

Abstract

Background

Tuberculosis is an air-borne contagious disease caused by slow-growing Mycobacterium tuberculosis (Mtb). According to Global Tuberculosis Report 2022, 1.6 million people were infected by tuberculosis in 2021. The continuing spread of drug-resistant tuberculosis (TB) is one of the most difficult challenges to control the tuberculosis. So new drug discovery is essential to the treatment of tuberculosis. This study aims to develop a QSAR model to predict the antitubercular activity of tetrahydronaphthalene derivatives. The QSARINS was used in this study to develop the QSAR predictive model.

Results

A number of tetrahydronaphthalene derivatives with MIC₉₀ values were obtained from the literature to develop the QSAR predictive model. The compounds were divided into two sets: a training set consisting of 39 compounds and a test set containing 13 compounds. The best predictive Model 4 has R²: 0.8303, Q²_LOO: 0.7642, LOF: 0.0550, Q²-F₁: 0.7190, Q²-F₂: 0.7067, Q²-F₃: 0.7938 and CCC_ext: 0.8720. Based on the developed QSAR equation, the new compounds were designed and subjected to molecular docking, molecular dynamics and ADMET analysis.

Conclusion

In the QSAR model, the molecular descriptors of MATS8s, Chi4, bcutv8, Petitjeant and fr_aniline were highly influenced the antitubercular activity. The developed QSAR model helps to predict the antitubercular activity of tetrahydronaphthalene derivatives.

1 Background

Mycobacterium tuberculosis is an air-borne bacterial infectious agent that causes tuberculosis (TB), a chronic disease that spreads through the air [1]. It is one of the most common causes of mortality from a single pathogenic agent that mainly affects the lungs. It is an acid-fast, slow-growing gram-positive bacteria with a very impermeable cell wall. According to Global Tuberculosis Report 2022, 1.6 million people were diagnosed with TB in 2021, an increase of 4.5% from 1.1 million people in 2020. In 2021, there were an estimated 14 lakh TB patient death among HIV-negative people and 1.87 lakhs TB patients’ deaths among HIV-positive people [2].

The World Health Organization (WHO) has issued guidelines for treating tuberculosis through a program known as Directly Observed Treatment (DOTS short-course). This approach involves a six-month treatment regimen of four first-line TB drugs: rifampicin, isoniazid, pyrazinamide and ethambutol. The treatment begins with a two-month course of all four drugs, followed by a four-month course of rifampicin and isoniazid [3]. One of the major issues with the treatment of TB is the emergence of multidrug-resistant tuberculosis. Genetic mutations mainly cause resistance [4]. The emergence of multidrug-resistant (MDR) and extensively drug-resistant (XDR), and totally drug-resistant-TB (TDR), along with the toxicity of existing antitubercular drugs leading to adverse side effects, leads to reviving the efforts of antitubercular drug discovery [5].

Therefore, it is imperative to develop new therapeutics with improved toxicity profiles in order to overcome these challenges. After more than 40 years, the discovery of delamanid (mycolic acid biosynthesis inhibitor) and bedaquiline (mitochondrial ATP synthase inhibitor) provided a degree of relief in the treatment of MDR-TB. However, the two medications have some notable side effects, such as ADME issues and hERG toxicity [6]. The discovery of new antitubercular agents with low toxicity is the only choice to face this difficult situation. But it is a costly and time-consuming process. Nowadays, various in silico tools are used to speed up the drug discovery process. So in this work, we used the quantitative structure–activity relationship (QSAR) to develop the predictive model. The QSAR is an adorable computational tool; it helps to find the relationship between the physiochemical properties of the compounds and its biological activity. Throughout the past thirty years, lead optimization has been sped up using the QSAR (quantitative structure–activity relationship) model. It helps to predict the pharmacological activity of novel compounds; it aids in the development of mathematical QSAR models. It gives ideas for lead optimization and virtual screening priority [7].

The current study deals with QSARINS, an adorable computational software that helps develop a robust QSAR model, which allows the design of novel compounds with higher mycobacterial tubercular inhibitory activity [8, 9]. In this study, we are investigating the in silico theoretical studies to develop novel mycobacterium tuberculosis inhibitors. We present here the results of a QSAR study on derivatives of tetrahydronaphthalene that are effective inhibitors of mycobacterium tuberculosis. Tetrahydronaphthalene is one of the excellent building blocks in medicinal chemistry. It has a wide range of biological properties, including antitubercular activity [10], anticancer [11, 12], anti-inflammatory [13], antiviral [14], anti-Alzheimer [15], antimicrobial [16] and antimalarial [17].

With the significant biological activity of tetrahydronaphthalene and its derivatives, we chosen this scaffold for our QSAR study. The current study aimed to construct a QSAR model of new tetrahydronaphthalene derivatives for their mycobacterium tuberculosis inhibitor and to explain how structural variation affects the antitubercular activity. The tetrahydronaphthalene derivatives showed good mycobacterium tuberculosis inhibition activity with a range of MIC₉₀ values of 0.21μg/mL to 13.5μg/mL. The QSAR model development was performed to investigate and explore the structural characteristics of tetrahydronaphthalene analogs required to inhibit Mycobacterium tuberculosis. According to the OECD regulatory guidelines, thorough validation was done, and its parameters were well-fit. The applicability domain was also developed. We performed ligand-based in silico studies to interpret the relationship between structural features and the biological activity of the compounds. In this paper, we present an analysis of 2D-QSAR studies conducted on the tetrahydronaphthalene derivatives that have been identified as effective inhibitors against mycobacterium tuberculosis [18,19,20].

2 Materials and methods

2.1 Softwares

QSARINS software was used to develop the 2D QSAR predictive model. The 2D structure of the compounds was sketched by ACD/Labs ChemSketch (Freeware) 2021.1.1. The energy minimization of the 2D structure was done using the MM2 force field in Chem3D pro version 12.0.2.1076 [21, 22]. PaDEL Descriptor software (version 2.20) was used to construct the molecular descriptors [23, 24]. The physiochemical and pharmacokinetic properties (ADME) of the compounds were evaluated using the SWISSADME online server. The quality of the selected protein was evaluated by the Ramachandran plot using a procheck webserver. The various toxicity studies including hepatotoxicity and carcinogenicity were predicted using the Protox-II webserver [39, 40]. The active site of the selected target protein was predicted by the Supercomputing Facility for Bioinformatics & Computational Biology (SCFBio), IIT- Delhi as DBT'S center of excellence. The molecular docking of the ATP synthase was performed using AutoDock Vina.

2.2 Molecular sketching and 3D optimization

QSARINS software was used to develop the MLR-QSAR predictive models by ordinary least squares. The development of QSAR was pioneered at the University of Insubria, Italy [25]. QSARINS uses a hybrid approach that combines genetic algorithms (GAs) and multiple linear regressions (MLRs) [25]. It helps build QSAR models that are both highly predictive and comprehensible [26]. The compounds structure was drawn by ChemSketch and saved as a mol file format. The energy minimization of the 2D structure was done by using the MM2 force field in Chem3D pro and saved as the mol file format. The molecular descriptor for 55 compounds was calculated by the PaDEL Descriptor version 2.20 software [27, 28].

2.3 Dataset collection, splitting and model development

For the studies, 55 molecules with antitubercular activity on H37Rv were collected from the literature. Wide variations in activity further enhance the performance of compounds with good activity [10]. The obtained data set was split into two sets, namely the training set or modeling set (70%) and the test set or prediction set (30%). The training set is used to model development and the test set is used to validation of the developed model. This work uses the MLR using the ordinary least-square approach to construct a good QSAR model. Multiple linear regression was used to correlate the relationship between the biological activity (MIC₉₀) and molecular descriptors. Here biological activity is the dependent variable and molecular descriptors are the independent variables. A genetic algorithm was used to selection of good molecular descriptors to build the QSAR model. It shows the direct correlation between independent variable X and dependent variable Y (descriptors). The dependent variable Y (biological activity) in an MLR analysis depends on independent variable X (molecular descriptors) [29].

$$Y \, = \, k1*x1 \, + \, k2*x2 \, + \, k3*x3 \, + \, C$$

Y is the dependent variable (biological activity), x is the independent variable (molecular descriptor), ‘C’ is a regression intercept and ‘k’s are the regression coefficients [30].

2.4 Applicability domain

The applicability domain of a QSAR model is a theoretical spatial region defined by the molecular properties or structural information of the compounds used in the QSAR model development. In which the William plot was used to define the applicability domain so that reliable prediction of the QSAR model has a leverage value (h) of the compounds below the critical leverage value (h*) with a standard deviation value of ± 3. The compounds that are outside the spatial region of the applicability domain may be considered an outlier. This method used to find the influential molecules reveals the model's ability to make predictions in the limited region. The term "outlier" refers to any molecule that deviates from the domain space or range. If the cross-validated standardized residual produced by the model is greater than ± 3, then the data are said to be an outlier. The leveraged method goes into more detail on the QSAR model's applicability domain [31, 32]. The formula for the warning leverage (H*) is as follows:

$$H* \, = \, 3\left( {p \, + \, 1} \right) \, / \, n$$

H* = The warning leverage value to check the effective molecule.

p = The number of molecular descriptors to build the QSAR model.

n = The number of molecules in the training dataset.

2.5 Y-randomization

It serves as an external validation factor to assess the validity and reliability of the developed QSAR model. The Y-randomization test was performed on the training dataset. In the Y-randomization test, the dependent variables (biological activity data) were shuffled, while the independent factors remained unchanged to create the MLR-QSAR model [33]. The developed QSAR model should have low R² and Q² values, in order to ensure that the model equation is reliable. R²_Yscr needs less to get a reliable model. It confirmed that the developed model was not obtained by the random chance correlation. If the developed QSAR model’s internal validation process was good, then the values of R² and Q² of each iteration, and their averages (R²_Yscr and Q²_Yscr) should be lesser with respect to the values of the model [34, 35].

2.6 Validation of the developed QSAR model

The stability, reliability, predictive strength, fitting criteria and the robustness of the QSAR models were validated by internal validation parameters and external validation parameters using QSARINS software. We used various statistical parameters like correlation coefficient (R²), adjusted correlation coefficient (R²_adj), leave-one-out squared correlation coefficient (Q²_LOO) and leave-many-out squared correlation coefficient (Q²_LMO) to evaluate the internal validation of the model. We also used some statistical parameters like Q²-F₁, Q²-F₂, Q²-F₃ and concordance correlation coefficient (CCC) to evaluate the internal validation of the model. Higher R² for goodness of fit, Higher Q²_LOO for robustness, lowest difference between R² and Q²_LOO for stability, low K_xx for less correlation among the descriptor and high K_xy for high correlation among the descriptor and biological response. The lower RMSE_Tr, RMSE_CV, and RMSE_ext must be less and as close as possible to get good model. The minimum suggested statistical values tabulated must be met for the QSAR model to be considered acceptable [36,37,38].

2.7 ADMET study

The physicochemical and pharmacokinetic properties (ADME) of the compounds were evaluated using SWISSADME online server. The druglikeness property was evaluated by Lipinski’s rule of five. According to the Lipinski rule of five, molecular weight should be less than 500 dalton, hydrogen bond donor should be less than 5, hydrogen bond acceptor should be less than 10 and partition coefficient (LogP) should be less than 5. The various toxicity studies including AMES toxicity, hepatotoxicity, hERG I inhibitor and hepatotoxicity were predicted using an online server called pkCSM server [39, 40].

2.8 Protein and ligand preparation

The drug target proteins were retrieved from the RCSB-PDB database with higher resolution. The crystal structure of the ATP synthase (PDB ID 7NK7) was predicted by the electron microscopy method, resolution of 2.11 Å [41, 42]. This protein is mainly involved in the catalysis for the synthesis of ATP from ADP and phosphate in mycobacterium tuberculosis. The retrieved proteins were prepared by removal of cocrystallized ligand and water molecules using Molegro Molecular Viewer and saved as the PDB file format. The 3D structure of ligand molecules was sketched by ChemSketch and saved in mol file format. These 2D structures of the ligand molecules were energy-minimized by Chem3D pro 12.0 software and saved as the PDB file format. The physiochemical properties of the molecules were calculated using SWISSADME online server [43,44,45].

2.9 Active site prediction and molecular docking

The active site of the selected target protein was predicted by the Supercomputing Facility for Bioinformatics & Computational Biology (SCFBio), IIT- Delhi as DBT'S center of excellence. The newly designed molecules were docked to the drug target protein. The molecular docking was to find the potential binding affinity of the designed compounds with the target protein. The process of molecular docking the ligands to the target protein encompasses several stages: preparation of the protein, preparation of the ligand, generation of the receptor grid, and actual molecular docking. The molecular docking study of the ATP synthase and tetrahydronaphthalene derivatives was performed using AutoDock Vina with grid coordinates of Center x (246.074), center y (194.992) and center z (251.111), respectively. The grid spacing was set to 70 Å with default docking parameters. The binding interaction of the protein and ligand was visualized and analyzed using Biovia Discovery Studio Visualizer 2021 [46, 47].

2.9.1 Molecular dynamic simulation

Molecular dynamics (MD) simulation was carried out using GROMACS 2022.2. The following steps were utilized.

2.9.2 Preparation of enzyme

The three-dimensional (3D) model of the ligand–protein complex was exported to.pdb format using Pymol. The dynamic behavior of the complex was evaluated using molecular dynamic (MD) simulation in the GROMACS package program (version 2022.2) [48,49,50]. Protein topology was constructed by pdb2gmx with the CHARMM27 force field [51], and ligand topology was generated using the SwissParam server [52].

2.9.3 Setting up a system for simulation

After applying the force field, the complex was inserted into the system. It was solvated with the TIP3P water model [53] in a cubic box greater than 1 nm from the edge of the protein with periodic boundary conditions. The system was neutralized by adding Na + ions, and energy minimization was done for 50,000 steps using the steepest descent algorithm. It was then followed by 100 ps of NVT simulation at 300 K and 100 ps of NPT simulation to equilibrate the whole system. Leapfrog algorithm was employed in the constant-temperature, constant-pressure (NPT) ensemble to separately couple each component like protein, ligand, water molecules, and ions [54]. The Berendsen temperature and pressure coupling constants were set to 1 and 2, respectively, to keep the system in a stable environment (300 K temperature and 1 bar pressure) [55]. The NPT ensemble step in molecular dynamics (MD) simulations plays a vital role by controlling pressure and temperature. This equilibrium-seeking step enhances the accuracy and stability of simulations, ensuring that they reflect the true thermodynamic behavior of molecules. It also facilitates efficient sampling of phase space, which is critical for obtaining statistically relevant results. Additionally, the NPT ensemble maintains realistic density and compressibility, which is crucial for systems with different phases or phase transitions. Finally, MD simulation for 100 ns was performed in isothermal and isobaric condition ensemble at 300 K. The pressure coupling with time constant was set at 1 ps to maintain pressure constant at 1 bar, and the LINCS algorithm [56] was used to constrain the bond lengths. The van der Waals and Coulomb interactions were truncated at 1.2 nm, and the PME algorithm [57] inbuilt in GROMACS was used to minimize the error from truncation.

2.9.4 Visualization and analysis of simulation

The trajectory file was visualized through VMD (Visual Molecular Dynamics) 1.9.2. [58] and analyzed by indigenously developed tool HeroMDAnalysis and Xmgrace 5.1.25 [59, 60].

2.9.5 Molecular mechanics generalized born and surface area (MM/GBSA) calculations

The binding free energies of the protein–ligand complex were calculated at the time of 0 ns, 50 ns and 100 ns by using molecular mechanics/generalized born surface area (MM/GBSA) methods. All the water molecules and ions were removed before MM/GBSA calculations. MM/GBSA free energy decomposition analysis was also used to calculate the binding free energy between the protein–ligand complex.

3 Results

3.1 QSAR model development and validation

The tetrahydronaphthalene and its antitubercular activity (pMIC₉₀) were used for this QSAR modeling study. The QSARINS software was utilized to generate multiple MLR-QSAR models. In this study, we focus exclusively on the most optimal QSAR models. For the studies, 55 molecules with antitubercular activity on H37Rv were collected from the literature. Wide variations in activity further enhance the performance of compounds with good activity. The molecular descriptors were calculated using PaDEL descriptor software and ChemDes online web software. It calculated around 2813 molecular descriptors for all the compounds. During the prereduction of the molecular descriptor step, it removed around 2413 molecular descriptors which have more than 80% of constant value and more than 95% of intercorrelated among the molecular descriptors. The compounds THN-13, THN-14 and THN-56 were entirely different from the rest of the compounds in the applicability domain. So they were excluded as outliers. The obtained data set was split into two sets, namely the training set or modeling set (70%) and the test set or prediction set (30%). The training dataset (39 compounds) was used to QSAR model development, and the test dataset (13 compounds) was used for validation of the developed QSAR model. Table 1 lists the various statistical parameters for developed QSAR models. Table 2 shows the correlation matrixes among the molecular descriptors in developed QSAR model 4. Table 3 enlists the calculated and predicted antitubercular activity by model equation and LOO method. The Williams plot helped to find the applicability domain of the developed QSAR model. The applicability domain was defined within a squared area of leverage threshold of h* = 0.423 (x-axis) and ± 3 standard residuals (y-axis) for QSAR model 4. h* = 3p´/n, where n is the number of molecules used in the developed QSAR model, p´ is the number of the molecular descriptors in the model plus one, and h* is the warning leverage value.

Table 1 Statistical parameters of QSAR model 1, 2, 3 and 4

Table 2 Correlation matrixes between the descriptors

Table 3 Calculated and predicted antitubercular activity by model equation and LOO

Figure 1 shows a scatter plot of the experimental value versus the predicted value by the model equation. In this diagram, the pale circles depict the compounds belongs to the training set, while the deep circles symbolize the compounds belonging to prediction set. Figure 2 depicts the scatter plot of experimental value versus predicted value by the leave-one-out (LOO) method. Figure 3 shows the scatter plot of leave-many-out (LMO) validation. This plot explains the intercorrelation among the molecular descriptors. Figure 4 illustrates the scatter plot of residual predicted by leave one out (LOO). Figure 5 shows the Williams plot of the Hat diagonal leverage value versus standardized residuals. It helps to find the outliers of the developed QSAR model. In this developed model, only one compound shown as an outlier. Figure 6 illustrates the scatter plot of the Y-scrambled model. This Y-scramble plot explains the external validation parameters.)

Fig. 1

Scatter plot of experimental pMIC₉₀ versus predicted value by model equation

Fig. 2

Scatter plot of experimental pMIC₉₀ versus predicted by leave one out (LOO)

Fig. 3

Scatter plot of Leave many out (LMO) validation compared with the original model

Fig. 4

Scatter plot of residual predicted by leave one out (LOO)

Fig. 5

Williams plot of Hat diagonal leverage value versus standardized residuals

Fig. 6

Scatter plot of Y-scrambled model compared with the original model

QSAR Model 1

$$\begin{aligned} {\text{pMIC}}_{{{9}0}} = & { 94}.{24}0{3 } + { 1}.0{546}*{\text{ MDEC}} - {11 } + { 2}.{9961}*{\text{ bcutv8 }} - \, 0.{8756}*{\text{ bcutm6 }}{-}{ 3}0.{7}0{53}*{\text{ bcute1 }} \\ + { 9}.{7121}*{\text{ bcutp1 }}{-}{ 17}.{48}00*{\text{ petitjeant }}{-} \, 0.{6645}*{\text{ fr}}\_{\text{aniline}} \\ \end{aligned}$$

$$\begin{gathered} {\text{R}}^{{2}} = \, 0.{91}0{3},{\text{ Q}}^{{2}}_{{{\text{LOO}}}} = \, 0.{8629},{\text{ Q}}^{{2}} - {\text{F}}_{{1}} : \, 0.{5776},{\text{ Q}}^{{2}} - {\text{F}}_{{2}} : \, 0.{559}0,{\text{ Q}}^{{2}} - {\text{F}}_{{3}} : \, 0.{69}00, \hfill \\ {\text{CCC}}_{{{\text{ext}}}} : \, 0.{8282} \hfill \\ \end{gathered}$$

QSAR Model 2

$$\begin{aligned} {\text{pMIC}}_{{{9}0 }} = & {94}.{6383 } + { 3}.{5}0{63}*{\text{ bcutv8 }} - \, 0.0{115}*{\text{ LogP2 }} - { 3}0.0{538}*{\text{ bcute1 }} + { 7}.{6727}*{\text{ bcutp1}} \\ \, - {15}.{267}0*{\text{ petitjeant }} - 0.{6416}*{\text{ fr}}\_{\text{aniline}} \\ \end{aligned}$$

$$\begin{gathered} {\text{R}}^{{2}} = \, 0.{8747},{\text{ Q}}^{{2}}_{{{\text{LOO}}}} = \, 0.{82}0{5},{\text{ Q}}^{{2}} - {\text{F}}_{{1}} : \, 0.{58}0{3},{\text{ Q}}^{{2}} - {\text{F}}_{{2}} : \, 0.{5618},{\text{ Q}}^{{2}} - {\text{F}}_{{3}} : \, 0.{692}0, \hfill \\ {\text{CCC}}_{{{\text{ext}}}} : \, 0.{81}0{6} \hfill \\ \end{gathered}$$

QSAR Model 3

$$\begin{aligned} {\text{pMIC}}_{{{9}0 }} = & {97}.{364}0 \, + \, 0.{1292}*{\text{ SHBint8 }} + { 2}.{3249}*{\text{bcutv8 }} - {3}0.{4}0{88}*{\text{ bcute1 }} + { 8}.{1841}*{\text{ bcutp1 }} \\ - {15}.{8477}*{\text{ petitjeant }} - 0.{623}0*{\text{ fr}}\_{\text{aniline}} \\ \end{aligned}$$

$$\begin{gathered} {\text{R}}^{{2}} = \, 0.{868}0,{\text{ Q}}^{{2}}_{{{\text{LOO}}}} = \, 0.{8274},{\text{ Q}}^{{2}} - {\text{F}}_{{1}} : \, 0.{58}0{4},{\text{ Q}}^{{2}} - {\text{F}}_{{2}} : \, 0.{5619},{\text{ Q}}^{{2}} - {\text{F}}_{{3}} : \, 0.{6921}, \hfill \\ {\text{ CCC}}_{{{\text{ext}}}} : \, 0.{71}0{7} \hfill \\ \end{gathered}$$

QSAR Model 4

$$\begin{aligned} {\text{pMIC}}_{{{9}0 }} = & {7}.{1879 } + { 1}.{1936}*{\text{ MATS8s }} + { 3}.{5538}*{\text{ bcutv8 }} - 0.{2939}*{\text{ Chi4 }} - {16}.{1}0{15}*{\text{ petitjeant }} \\ - 0.{4963}*{\text{ fr}}\_{\text{aniline}} \\ \end{aligned}$$

$${\text{R}}^{{2}} : \, 0.{83}0{3},{\text{ Q}}^{{2}}_{{{\text{LOO}}}} : \, 0.{7642},{\text{ Q}}^{{2}} - {\text{F}}_{{1}} : \, 0.{719}0,{\text{ Q}}^{{2}} - {\text{F}}_{{2}} : \, 0.{7}0{67},{\text{ Q}}^{{2}} - {\text{F}}_{{3}} : \, 0.{7938},{\text{ CCC}}_{{{\text{ext}}}} : \, 0.{872}0$$

3.2 Design, molecular docking and ADMET results

Based on the QSAR model, we designed some structural analogs of the tetrahydronaphthalene compounds, and its antitubercular activity was predicted. The Mycobacterium smegmatis ATP synthase F1 state 1 (PDB ID:7NK7) was chosen as the drug target, which is determined by the electron microscopy method with a resolution of 2.11 Å. The quality of the protein was determined by Ramachandran plot using the Prockeck webserver. Figure 7 depicts the Ramachandran plot of the ATP synthase protein. Based on the developed QSAR model, we designed 10 compounds enlisted in Table 4. All the newly designed compounds showed good docking scores ranging from − 7.5 kcal/mol to − 9.1 kcal/mol. The standard drug Bedaquiline had a docking score of − 7.5 kcal/mol toward the protein. The THN_RD9 had docking score of − 9.1 kcal/mole also had various Hydrogen bond interactions with ARG A:174, LYS A:175, LYS A:178 and GLN A:435, Pi-Cation interaction with LYS A:276 and Pi-Alkyl interaction with ALA A:180, VAL A:183, ALA A:222 and VAL A:219 amino acids of ATP synthase protein. Figure 8 illustrates the 3D and 2D interaction of THN_RD9 with ATP Synthase protein. All the designed compounds don't have any Lipinski violation and toxicity (hepatotoxicity and carcinogenicity) (Tables 5, 6, 7).

Fig. 7

Ramachandran plot of the ATP synthase

Table 4 The structures and SMILES of newly designed compounds

Fig. 8

3D interaction of THN_RD9 with ATP synthase protein

Table 5 ADMET properties of newly designed compounds

Table 6 Molecular docking result of the designed compounds

Table 7 Binding free energy of protein–ligand complex using MM/GBSA calculation

3.3 Molecular dynamic simulations results

In order to understand the conformational changes and evaluate the binding of THN_RD9 against Mtb ATP synthase (PDB ID: 7NK7), we have carried out MD simulation for a period of 100 ns for THN_RD9-Mtb ATP synthase complex. The simulation was evaluated using various statistical parameters including root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), h-bond interactions, and its % occupancies over the time. Figure 9 illustrates the protein–ligand complex of Mtb ATP synthase protein with THN_RD9 for molecular dynamics simulation.

Fig. 9

Protein–ligand complex of THN_RD9-Mtb ATP synthase, where protein is shown in cartoon representation and the ligand is shown in CPK representation

3.3.1 RMSD analysis

Analyzing the RMSD can give insights into any structural conformation that protein and ligand undergo during the simulation. The multiplot for protein Cα versus time for the simulation is shown in Fig. 10. Both the ligand and protein in the complex form have attained a plateau in RMSD values (of less than 0.3 nm), which indicated that the ligand–protein complex was stable during the simulation.

Fig. 10

Graphical representation of the plot showing protein Cα and Ligand RMSD (nm) versus time (100 ns) for (A) Mtb ATP synthase protein (green in color), and (B) THN_RD9 ligand (blue in color)

3.3.2 RMSF analysis

The protein-RMSF is useful for characterizing local changes along the protein chain. The plot for protein-RMSF (nm) versus residue number index is shown in Fig. 11. Notably, the plot describes a fluctuation of less than 0.3 nm for the protein. Particularly the residues, which were involved in binding with the ligand, have shown very small fluctuations. It also indicated the stability of ligand–protein in complex form.

Fig. 11

Graphical representation of the plot showing the protein RMSF (nm) versus residue index number of protein for THN_RD9-Mtb ATP synthase complex

3.3.3 H-bond interaction

Molecular interactions, particularly the h-bond interactions, are distance and angle depend and liable to disrupt under dynamic conditions. Herein, we have analyzed the ligand–protein interactions. The plot for the number of hydrogen vs time is shown in Fig. 12. From the plot, it was observed that THN_RD9 displayed 1–2 h-bond contact in binding Mtb ATP synthase. To access the residues involved in such interactions and their stabilities, the %occupancies vs the residues were also calculated.

Fig. 12

Pictorial representation of the number of h-bond contacts formed by THN_RD9 in complex with Mtb ATP synthase (PDB ID: 7NK7)

Figure 13 represents the histogram of % occupancies of the h-bond contacts formed by THN_RD9 in binding Mtb ATP synthase. This graph has displayed the ability of THN_RD9 to form stable h-bond contacts with residues VAL139, LYS212, LYS175 and THR215, with occupancies of 16.16, 6.86, 1.59 and 1.04%, respectively. So overall, it can be concluded that the THN_RD9 can be efficient in binding with Mtb ATP synthase.

Fig. 13

Histogram representation of %occupancies of the h-bond protein ligand contacts of THN_RD9 in complex with Mtb ATP synthase (PDB ID: 7NK7)

3.4 MM/GBSA calculation

Further, to understand the energy contributions of individual residues inside the active site, the per-residue interaction energies were calculated for the Mtb ATP synthase in complex with THN_RD9. The major contributing residues include ASP241, ASP272 and GLU331 (Fig. 14) The binding free energy of the protein–ligand complex was calculated using MM/GBSA calculation. In this result, the binding free energy of the initial time (0 ns) of the protein-THN RD9 complex was − 31.57 kcal/mol, at 50 ns the binding free energy was -39.85 kcal/mol; at the final stage(100 ns), the binding free energy was -34.29 kcal/mol.

Fig.14

Per-residue energy contribution using molecular mechanics/generalized born and surface area continuum solvation (MM/GBSA)

4 Discussion

We created four QSAR models that exhibited strong model fitting and adhered to the OECD guidelines for producing reliable predictive QSAR models. However, we selected model 4 as the best model based on the external validation parameters. Hence, we are focusing on the validation, significance and predictive ability of the developed QSAR model 4. This model exhibits strong internal predictive capability, as indicated by a Q² value of 0.7642. The goodness of fit of the developed model was the correlation coefficient (R² = 0.8303). The developed QSAR model’s predictive ability was evaluated by an external correlation coefficient (R²_adj = 0.7751). The lower R²_Yscr values indicate that the proposed QSAR models are not obtained by the chance correlation. The importance of the quantity of molecular descriptors in the developed QSAR model was validated by the minimal difference of 0.0257 (less than 0.3) observed between R² and R²_adj. The stability of the developed QSAR model was confirmed by the lowest difference of 0.0661(less than 0.3) between the R² and Q²_LOO. In this method, we excluded one compound and computed the model with the remaining compounds and then we predicted the activity of the excluded one. The next stronger technique included in the QSARINS is the leave-many-out (LMO) method, in which 30% of the compounds were left out to develop the predictive model and to study the behavior of the model.

All the statistical parameters of the four models are given in Table 1, in which the fourth model’s coefficient of determination, denoted as R², is 0.8303, which helps to evaluate the goodness of fit. A good model should have an R² value greater than 0.6. In the fourth model, the coefficient of determination denoted as R² is 0.8303, aiding in the assessment of goodness of fit. A robust model should exhibit an R² value greater than 0.6. R²_adj represents the adjusted coefficient of determination; it gives insight into the suitability by addition of new molecular descriptors to the QSAR model. For a model to be considered good, R²_adj should exceed 0.6. LOF refers to Friedman’s lack of fit criteria, which assesses overfitting in the QSAR model. For a model to be deemed good, LOF should be below 0.3. Kxx represents the overall correlation among the descriptors. Delta K indicates the difference in correlation between the descriptor (kx) and the correlation difference between the descriptor and the response (kxy).

RMSE_tr stands for root-mean-square error in the training set calculations. For a model to be considered good, its RMSE should be under 0.3. MAE represents the mean absolute error in the adjustment calculated within the training series; for a model to be considered good, it should be below 0.3. S stands for the standard estimation error. A good model should exhibit an S value lower than 0.3. Additionally, both R²_Yscr and Q²_Yscr should be minimized to indicate that the developed model is not obtained by chance correlation. The model 4 fulfills all the required statistical parameters. Q²-F₁, Q²-F₂, and Q²-F₃ values are helpful to quantify the predictive power of the developed model; it should also exceed 0.7 for a good predictive model. The CCC represents the concordance correlation coefficient; for a model to be considered good, it must surpass 0.85. This CCC value quantifies the similarity of the predicted and experimental values. While Q²F₃ and CCC are the most stable and reliable parameters. The value of F denotes Fisher’s F value, which should exceed the theoretical threshold. The domain of applicability of the selected model was illustrated in William's plot. It showed that most of the compounds were placed within the theoretical space of the model's applicability domain. One compound (THN-42) with leverage values (h = 0.6130) is greater than the warning leverage value (h* = 0.462), However, the plot shows that standardized residual value was within the limits. This compound must be considered as possible outliers; however, its behavior is relatively similar to the remaining compounds, so we did not exclude them from this study.

The Organization for Economic Co-operation and Development (OECD) Community gave five important principles to develop a good QSAR predictive model. OECD principle 1 deals with a defined end point which means it refers to a pharmacological or biological activity that could be measured. OECD principle 2 associated with unambiguous algorithm, it explained about how the molecular descriptor of the compounds there biological activity are related. OECD principle 3 deals with the Applicability Domain. It aids in identifying outliers, which refer to compounds that fall outside the spatial region of the applicability domain. OECD principle 4 deals with robustness, predictive ability and goodness of fit of the developed QSAR model. OECD principle 4 deals with the mechanistic interpretation, which is not mandatory for the model.

In model 4, various molecular descriptors such as MATS8s, bcutv8, Chi4, Petitjeant and fr_aniline influence the antitubercular activity, in which MATS8s and Chi4 were positively correlating and contributing toward the antitubercular activity. Other descriptors like bcutv8, Petitjeant and fr_aniline were negatively contributed to the antitubercular activity, in which bcutv8 was related to the atomic van der Waals volume. The molecular descriptor Chi is the atomic connectivity indices. MATS8s descriptor is the Moran autocorrelation of lag 8 weighted by I-state. Petitjeant is the topological descriptor. Based on the developed QSAR model equation, the newly designed 10 tetrahydronaphthalene compounds were subjected to molecular docking and ADMET analysis. All the designed compounds had good docking scores (greater than -7.1 kcal/mol) and zero Lipinski violation. The THN-RD09 compound has high pMIC_90, good binding affinity, zero Lipinski violation and no toxicity.

5 Conclusion

In the present study, the QSAR study was carried out on a series of tetrahydronaphthalene derivatives with antitubercular activity. The molecular descriptors were calculated using by the PaDEL Descriptor software and ChemDes webserver. It helps to develop a stable, predictive and robust model with fulfilled OECD principles. The developed QSAR model's robustness was evaluated through various internal (Q²_LOO, Q²_LMO and R²_Yscr) validation parameters, and its predictive ability was assessed using external validation. In the QSAR equation, the molecular descriptors of MATS8s and Chi4 positively correlated with antitubercular activity. While the bcutv8, Petitjeant and fr_aniline were negatively correlated with antitubercular activity. Modification of these molecular descriptor parameters will lead to good antitubercular activity. The Mtb ATP synthase inhibition was greatly dependent on atomic van der Waals volume, Moran autocorrelation, connectivity indices, topological factor and a fragment of aniline. Based on the developed QSAR model equation, the newly designed 10 tetrahydronaphthalene compounds were subjected to molecular docking and ADMET analysis. This developed QSAR predictive model will help to predict the antitubercular activity of tetrahydronaphthalene derivatives. The molecular dynamics (MD) simulations were performed for 100 ns, using the Gromacs package. It evaluated the conformational stability and alteration of protein–ligand complexes during the simulation. Thus, our findings confirmed that newly designed THN-RD9, a tetrahydronaphthalene derivative and known strong ATP synthase inhibitor, may be used as an important “lead” molecule to be developed as an antitubercular drug in the future. However, we need to do an in vitro and in vivo study to confirm its antitubercular activity.