Machine learning approaches in designing anti-HIV nitroimidazoles: 2D/3D QSAR, kNN-MFA, docking, dynamics, PCA analysis and MMGBSA studies

2.1 Dataset and biological activity

The QSAR analysis was conducted using a series of 20 compounds that were substituted for nitroimidazole as highly selective and potent inhibitors of non-nucleoside reverse transcriptase (Alderton et al., 2001). To associated with the structural characteristics (physicochemical descriptors), the compound’s logEC₅₀ (pEC₅₀) values were converted from their micro molar growth inhibitory concentrations. Table 1 displayed the chemical structures of 20 such compounds together with estimates of their biological activity used for generation of 2D and 3D QSAR models.

Close

Export to PPT

2.2 Dataset division into training and test sets

Compounds were sketched out and transformed into 3D structures using the 2D draw program. The Merck molecular force field and atomic charges were used in the batch energy minimization approach to do energy minimization and geometry optimization. The medium's dielectric constant was 1, the maximum number of cycles was 1000 and the convergence criterion (RMS gradient) was 0.01. To ensure the accuracy of our QSAR model we conducted a search on each energy minimized structure. This involved employing a strategy to the RIPS approach, where we introduced perturbations to the coordinates of each atom in the molecule generating fresh molecular conformations. When creating the training and test sets for our model we took into account the importance of having structures and a range of biological activities similar to those in the training set. To achieve this we utilized the sphere exclusion method (Bhadoriya et al., 2013; Bhatia et al., 2017) which determined the radius based on dissimilarity values (as shown in Table 1). This method allowed us to divide our dataset into two groups; a training set and a test set. Each group consisted of 20 molecules. Furthermore, we randomly selected 20 chemicals for our training set and 6 compounds for our test set. This division was essential in constructing and validating our model.

2.3 Molecular modeling

2.4 3D-QSAR analysis

2.4.5

2.4.5 Alignment rules

The molecular alignment of the supplied set of molecules was done to observe the structural diversity within them. A typical rectangular grid was then created enclosing the molecules. Unsubstituted nitroimidazole (the template structure) was employed for alignment while taking into account the shared parts of the series, as illustrated in Fig. 2A and 2B. The reference molecule as shown in Fig. 2A was chosen as a template molecule because of its potent inhibitory activity which made it a reliable lead molecule. The template alignment approach was used to superimpose every molecule from the series following optimization using the template structure and the reference molecule. The template-based alignment was applied to appropriately align the given set of molecules in accordance with the kNN-MFA methodology after optimization.

Close

Fig. 2

(A) Template molecule used for QSAR study and (B) Stereoview of all 20 aligned molecules.

Export to PPT

2.5 Molecular docking studies

The co-crystallized protein structure was obtained from the publicly accessible protein data bank (pdb id: 1HIV) (http://www.rcsb.org) whereas molecular docking investigations were carried out using a chosen ligand and reference molecule (Deeks et al., 2015). Then, using the protein preparation wizard the downloaded protein was subjected to different steps of adding the missing H atom and correctly allocating the bond order. Protein co-factors received formal charges, and an internal ligand heterostat was created. The last step was to perform OPLS2005 forcefield energy minimization to eliminate any unphysical interactions between the protein or ligand atoms and maintain the relaxed state of the co-crystallized protein complex. Any water molecules within 5 were also removed. The receptor grid was created in order to pinpoint the precise location where the compounds under study should bind. The centroid of the co-crystallized crystals was located within a 10 Å site, and default parameters were used with a 1.0 scaling factor. The designed nitroimidazole ligands with the reference molecule were put through the ligprep tool of Maestro utility. With the use of the OPLS2005 forcefield, the potential conformers, tautomers or stereoisomers of the ligands were created. For molecular docking investigations using the XP mode and a Van der scaling of 0.8, these synthesized ligands with energy minimization were employed. In the resulting grid's centroid spaces, 5 poses for each ligand were then captured and the best scoring pose was examined using XP visualizer.

2.6 MD simulation studies

Using the Desmond V3 program created by the DE Shaw group, chosen compounds with a satisfactory docking score were simulated. Three step procedure using the system builder, minimization and molecular dynamics run were employed to produce the trajectory file. The trajectory file was subsequently examined to get the required outcomes. The integrated ligand–protein complex discovered through docking experiments was single point charged (SPC) and a solvated orthorhombic border box with a distance of 10 Å was created. The produced system was neutralized using Na⁺ and Cl⁻ ions and a physiological salt content of 0.15 M was maintained. The system was then put through an energy minimization process utilizing NPT at 300 K and 1 bar of pressure. A 100 ns MD simulation was performed on this balanced system at standard temperature and pressure of 1 bar in order to produce the trajectory file. The simulation interaction diagram utility of the Desmond module was used to import the .cms file after a successful MD run and the outcomes were examined. We created RMSD (root mean square deviation), RMSF and ligand 2D interactions plots in order to examine the stability of complexes.

2.7 Principal component analysis (PCA) analysis

Principal component analysis (PCA) was employed to get a deeper understanding of the ways in which ligand binding impacts the Cα atom dynamics in ttbk1. Using this statistical method, we were able to understand the relationship between the different Cα conformations sampled throughout the four simulated systems throughout the molecular dynamics (MD) simulations and explain the Cα root-mean-square deviation (RMSD). The Bio3D R software tool developed especially for comparing protein structures made the study easier. Essentially, we used VMD to import the Desmond trajectory files (.dtr) and CHARMM/NAMD trajectory format (.dcd). files were made. We first overlaid the Cα atoms from every trajectory frame onto the minimised structure PDB file in order to get the data ready for PCA. After completing this procedure, we were able to obtain numerical coordinates (x, y, z) of the Cα atoms, which were imported into RStudio along with the relevant.dcd files. A lower-dimensional representation of the structural dataset was created by projecting the minimised structure and MD trajectory snapshots onto the subspace defined by the largest principal component (PC), which represents the highest variance in Cα atom locations. By using the pca.xyz function in the Bio3D package, we were able to project this image and capture all most important conformational changes that we saw throughout the simulations.

2.8 Binding free energies calculations (MMGB/SA calculations)

The computational technique known as MMGB/SA (Molecular Mechanics/Generalized Born Surface Area) is used in computational chemistry and structural biology to determine the binding free energy between proteins and ligands (small molecules). In drug development and studies of protein–ligand interactions, it is commonly used to predict the binding affinity of potential drug candidates to target proteins. The generalized Born solvent models, surface area calculations and molecular mechanics which makes use of force fields to describe a system internal energy are all combined in the MMGBSA method. The following is an equation (1) of the binding free energy (G bind):

2.9 In-silico ADME calculation

Many Drugs Design & Discovery project have been failed because of inadequate pharmacokinetic profiling of many potentially effective compounds. Even few molecules have been withdrawn from the market because of severe toxicities during its clinical usage leading to time loss and financial burden on any organization. Therefore, early-stage calculation ADME properties is important for successful drug development. In present study, QikProp (version 3.0) from Schrödinger, LLC was used for predictions virtual ADME. The parameters like drug-likeness of the ligand molecules based on Lipinski's Rule of Five and Rule of 3, lipophilicity, blood–brain barrier permeability, molecular weight and the numbers of H bond donors and acceptors etc.

3.1 2D QSAR (2D QSAR for III_B strain)

2D QSAR studies play a crucial role in drug design project by correlating the structural features of chemical compounds with their biological activities. These studies mainly focus on different structural descriptor such as hydrophobicity, electronic properties, atomic correlation and steric factors with their biological activities. These values can be easily calculated by studying about the 2D QSAR analysis and these models are less expensive, rapid and robust making them an essential tool of early-stage drug design or lead optimization (Chtita et al., 2022; Chtita et al., 2022). Here, the atom based 2D QSAR study of reported nitroimidazole was performed and results were discussed.

3.1.1

3.1.1 Interpretation

According to the equation, our model explains a significant 92.41 % (r² = 0.9241) of the total variance in the training set, according to the established equation (Table 2). With q² and pred_r² values of 77.06 % and 40.80 %, respectively, it displayed the excellent internal and external prediction ability. These results support the capability of generated models to give predictions that can be used to a wide range of situations. The QSAR equation's coefficients provide important information about the molecules that affect the biological activity. It is noteworthy that the negative coefficient for T_2_Cl_7 reveals the number of double-bonded atoms in a molecule that are distant from the chlorine atom by seven bonds, indicating a decrease in inhibitory activity for such configurations. On the other hand, a positive coefficient for SssSEindex, which stands for Electrotopological State Indices (ESI), shows that lower values of this descriptor are associated with stronger inhibitory action. The function of descriptors like T_N_Cl_6 (N atom bonded by single or double bond, separated from any oxygen atom) and SaasCE-index (number of C atoms with one single bond and two aromatic bonds) in anticipating inhibitory activity is revealed by further study. Increased inhibitory activity is shown by a larger T_N_Cl_6 value, whereas decreased inhibitory effects are indicated by a greater SaasCE-index value. It has been interesting to incorporate Electrotopological State Indices into our model for the prediction of inhibitory activity as shown in Fig. 3. These results lay a platform for understanding the structural components needed for biological activity, enabling rational drug design and development.

Table 2 Statistical Data obtained after the analysis of 2D-QSAR models of nitroimidazole derivatives.

Parameters	r2	q2	r2se	q2se	Pred_r2	Pred_r2se	F test
1(Model-1)	0.9241	0.7706	0.2017	0.3506	0.4080	0.4891	27.38
2(Model-2)	0.9528	0.8592	0.1377	0.2378	0.3358	0.6647	45.40
3(Model-3)	0.9047	0.6959	0.2353	0.4204	0.6540	0.3544	21.35

Mode-l (Test set: 14, 18, 19, 2, 20, 6).

pIC₅₀(column) = −0.7398(T_2_Cl_7) 0.7873(SssSEindex) −0.8244(T_N_Cl_6) 0.1539 (SaasCEindex) 1.2286.

Statistics: [n = 14; Degree of freedom = 9; q² = 0.7706; r²se = 0.2017; F test = 27.38; r² = 0.9241; q²se = 0.3506; pred_r²se = 0.4891; pred_r² = 0.4080].

Close

Fig. 3

(A) The 2D QSAR contribution plot shows us about the different QSAR parameter by using them we can design new potent derivatives on the basis of rationality, (B) data fitness plot shows the actual vs predicted activity of the best 2D QSAR model, and (C & D) Training and test set shows how the activity fluctuate from center point from actual to predicted, (2D QSAR using IIIrd B strain).

Export to PPT

3.2 2D QSAR for ROD strain

3.2.1

3.2.1 Interpretation

The outstanding external (pred_r²) and internal (q²) prediction abilities of the produced models which stand at 48.94 % and 82.99 %, respectively, were used to demonstrate their predictive powers. A noteworthy capacity to explain variance within the training set is also demonstrated by the model which accounts for a significant 94.12 % (r² = 0.9412) of the overall variance. The F-test result with a value of 36.027 is higher than the predicted threshold supporting the model's statistical significance (Table 3). SaaNEindex has a negative coefficient value [Electrotopological state indices for the number of N atoms bonded with two aromatic bonds]. Positive coefficient value of T_C_N_4 for C atoms which are four bonds away from N atom, either single or double bonded] indicated the lower values would give the potent inhibitory activity. A negative coefficient of T_T_N_5 results in a reduction in inhibitory activity, whereas a positive coefficient enhances inhibitory activity, according to the biological activity. A T_T_N_5 here indicate that single or double bonded N atoms which are five bonds away from another N atom in a molecule. A lower or positive coefficient value for T_2_Cl_7 indicates a reduced inhibitory effect on the molecule, whereas a negative coefficient value indicates increased inhibitory activity. T_2_Cl_7 term explained that any double bonded atom separated by 7 bonds from the chlorine atom in a molecule. The results analysis suggested that greater biological activity results in improved inhibitory activity, while lower biological activity results in decreased inhibitory activity. The obtained results reveal important molecular characteristics that affect biological activity. These findings suggested the important insights for design and discovery of newer molecule by indicating that increased biological activity was correlated with improved inhibitory activity and lower biological activity was related with decreased inhibitory activity as shown in Fig. 4.

Table 3 Statistical evaluation of 2D-QSAR models of nitroimidazole derivatives.

Trials	r2	q2	r2se	q2se	Pred_r2	Pred_r2se	F test
1(Model-1)	0.9412	0.8299	0.1698	0.2889	0.4894	0.5197	36.027
2(Model-2)	0.9649	0.8406	0.1249	0.2660	0.3883	0.6580	61.76
3(Model-3)	0.9464	0.9113	0.1350	0.1736	0.3919	0.7489	39.70

Model-(Test set: 1, 14, 17, 5, 8, 9).

pIC₅₀(column) = -12.6359(SaaNEindex) + 0.2926(T_C_N_4) −0.1239(T_C_N_5) −0.4142(T_2_Cl_7) 54.9262.

Statistics: [n = 14; r² = 0.9412; Degree of freedom = 9; q² = 0.8299; r²se = 0.1698; F test = 36.027; pred_r² = 0.4894; q²se = 0.2889; pred_r²se = 0.5197].

Close

Fig. 4

(A) The 2D QSAR contribution plot shows us about the different QSAR parameter, by using them we can design new potent derivatives on the basis of rationality, (B) data fitness plot shows the actual vs predicted activity of the best 2D QSAR model, and (C and D) Training and test set shows how the activity fluctuates from centre point from actual to predicted.

Export to PPT

3.3 3D QSAR (3D QSAR for III_B strain)

3D QSAR studies are the extension of 2D QSAR by adding the spatial arrangement of atoms and its correlation with biological activities. It gives the more comprehensive understanding how molecular conformations affect their biological activity. This allows for the lead optimization or molecular interactions optimization between drug and receptor leading to more potent and selective drug candidate. Considering these facts here, 3D QSAR study of reported nitroimidazole was done and results were reported Table 4. The values of k (3), q2 (0.7183), pred_r2 (0.7585), q2_se (0.3277), and pred_r2 se (0.3040) in Model 1 demonstrated the statistical significance of the QSAR equation thus derived and demonstrate the model predictive power of 75.85 % (internal validation) and 71.83 % (external validation) (Table 4). Moreover, the obtained values for Electrostatic field, −E_1656 (−0.8562, −0.6964) suggested that less bulky substituent groups are preferable in that location due to the fact that an increase in activity is favorably influenced by negative electrostatic potential. Similarly, Steric field value −S_1383(−0.045, −0.0866) indicated that a decrease in the activity is encouraged by a negative Steric potential, hence less bulky substituent groups are favored there. In this work, we used a thorough methodology to evaluate the predictive capability of our QSAR model. The projected inhibitory activity for the training and test sets are clearly displayed in Table 1. With the help of this division, we were able to evaluate how well our model was able to predict outcomes for new data. The data fitness plot for Model 1 is also shown in Fig. 5, giving a visual picture of the model's performance. A helpful view on the model's efficacy during training and its capacity to extend predictions to an external test set is provided by the observed vs predicted activity diagram shown in Fig. 5. The accuracy with which our model can predict inhibitory actions for compounds seen during training as well as those encountered for the first time in the external test set confirms the robustness of our approach. This demonstrates the usefulness and dependability of our QSAR model in applications for drug discovery and design.

Close

Fig. 5

(A) 3D alignment of molecules with the steric, (B) point data fitness plot shows the actual vs predicted activity of the best model, and (C and D) Training and test set shows how the activity fluctuate from center point from actual to predicted.

Export to PPT

3.4 3D QSAR for ROD strain

Owing to the importance of 3D QSAR for optimization of molecular interactions between ligand and receptor, the 3D QSAR studies of reported nitroimidazole derivatives for ROD Strain were determined and results were reported in Table 5. The statistical analysis of Model 1 gave its robustness and predictive power. The key metric highlights the importance of the derived QSAR equation, including k (2), q2 (0.9039), pred_r² (0.5622), q²_se (0.2071), and pred_r²_se (0.3355). The values obtained for Electrostatic field, −E_1812(−0.274, −0.265) shown that a less bulky substituent group is favored in that location and that a negative electrostatic potential is favorable for an increase in activity. It has been predicted that a negative electrostatic potential is desirable for an increase in activity, hence lower bulky substituents are preferred there. Similarly Steric field, −S_1120 (−0.1017, −0.0875); since it is suggested that a negative electrostatic potential is advantageous for an increase in activity, there is a desire for a lower substituent group there. The models have the good capacity to explain the variance in the training dataset accounting for almost 90.39 % of the observed variability, was explained by the high q² score of 0.9039 in internal validation. Moreover, the remarkable predictive performance of models on unobserved data is highlighted by the pred_r² score of 0.5622 in external validation, which explained that it can anticipate about 56.22 % of the variability in inhibitory activities of test set. The predicted inhibitory activity for both the training and test sets are summarized in Table 6 and displayed in Fig. 6 to give a clear insight of how well our model performed. The observed vs. anticipated activity diagram, confirms that our model is successful in identifying the underlying trends in the data. Due to its capacity to generalize predictions and direct decision-making processes, Model 1 is a significant asset in drug discovery and design after being subjected to this complete review.

Close

Fig. 6

(A) 3D alignment of molecules with the steric, (B) point data fitness plot shows the actual vs predicted activity of the best model, and (C and D) Training and test set shows how the activity fluctuates from center point from actual to predicted.

Export to PPT

3.5 Molecular docking studies

The molecules with good predicted IC₅₀ values as displayed in Table 7 were used further confirm the proposed hypothesis by means of computational study. Before the start of docking study, the software was validated by extracting and redocking of internal ligand without change in its heterostate. Then RMSD value compared to its original position was determined. The lower RMSD value 1.533 Å indicated that the software is validated and it is ready to use for further computational study. To predict how molecules, such as small novel HITs or ligands, interact with a target macromolecule, usually a protein, researcher utilize a computational technique called molecular docking. It is predicting the binding mechanism and binding affinity of a ligand to a protein receptor. This knowledge is essential for comprehending molecular interactions, creating novel medicines and improving ones that already exist. The downloaded protein 1HIV was a co-crystallized structure of HIV-1 protease. The structural analysis of this co-crystallized protein suggested that the protein known as Kinemage 2, is made up of two 99-residue polypeptide chains. While the second chain's numbers range from 101 to 199, the first chain's numbers range from 1 to 99. U-75875, a diol inhibitor with the numbers 201–206, is an additional compound. This inhibitor is situated within a well-defined and large active site cleft, just like other protease-inhibitor complexes described by Wlodawer et al. in 1992. Two separate “flaps” originating from each monomer make up the 23 Å long cleft that runs across the dimer interface on one side of the molecule. The HIV protease (1HIV protein) is an important target for the design and development of antiretroviral drugs (Thanki et al., 1992). Hence, we considered the molecular docking to examine the binding interactions between the identified HITs & this protein the results were compared the reference molecule. It was done to determine the binding affinities and mechanisms of action of identified with the active site of the 1HIV protein. The interactions between the identified HITs and the 1HIV protein were assessed by the molecular docking studies, which highlighted the importance of certain hydrogen bonds, hydrophobic contacts and electrostatic interactions. The results obtained were summarized in Table 8 and 3D and 2D interactive diagrams of reference molecule and identified HITs were displayed in Figs. 7 and 8, respectively. The docking scores of reference molecule and ZH-1, ZH-2 and ZH-3 were found to be −4.399, −4.712, −2.172 and −1.454 respectively. The overall good docking scores displayed by the ZH-3 might be due to the two H bond with Asp29 residue. The hydrophobic interactions were also observed with Pro81, Val82, Ile84, Leu23, Ala28, Ile50, Ile47, Val32, Ala28 and Ile84 active site residues with reference molecule and identified HITs. But in case of compounds ZH-3 additional hydrophobic residue Phe83 were observed. Similarly, the electrostatic −ve charged interactions were observed with Asp25, Asp29 and Asp30 residues in all the compounds including reference also. Two additional interactions were observed in case of ZH-3, firstly Pi-cation interactions with Arg8 and polar interaction with Thr80. Overall, these interactions with active site residues were quite identical in case of reference molecule and identified HITs making them good candidate for design of anti-retroviral therapy. Moreover, additional interactions observed in case of ZH-3, made this compound better candidate as development of anti-HIV drugs. However, these results were only the part of preliminary screening between receptor and ligands, the validation of the obtained docking results were needed with advance computational tools. Hence molecular dynamics (MD) simulations studies were performed to check the stability of docked complexes and observed interactions at atomic level.

Table 7 Newly Designed molecules for IIIB and ROD cell line with their predicted activities using the generated 2D and 3D QSAR models.

Sr. no	New molecules	Activity IIIrd	ROD
1		0.0004	0.103
2		0.050	0.08
3		0.084	0.0014
4		0.083	2.08
5		0.0169	0.096
6		0.0102	1.62
7		0.016	0.479
8		0.010	0.00098
9		0.0072	0.0030
10		0.00020	0.00427

Close

Fig. 7

3D interactive diagrams of different molecules under the study at the active site of 1HIV (A) reference molecules, (B) ZH-1, (C) ZH-2 and (D) ZH-3.

Export to PPT

Close

Fig. 8

2D interactive diagrams of (A) reference molecules, (B) ZH-1, (C) ZH-2 and (D) ZH-3 at the active pocket of HIV protease enzyme.

Export to PPT

3.6 MD simulations studies

All biological systems are thought of as dynamic networks of molecular interactions, in contrast to the singular rigid picture that molecular docking provides of the interaction between proteins and ligands (Bunker and Róg, 2020; Gelpi et al., 2015). We used the SPC model to simulate the dynamical behavior of the three promising docked complexes for 100 ns each in order to investigate the different conformations that these complexes might acquire in the solvated phase. For the 100 ns simulation run produces the 1000 frames (trajectories) of the protein–ligand complex structures, the protein C-atoms RMSD, the ligand RMSD with respect to protein, the Root Mean Square fluctuation (RMSF) and the protein ligand contact analysis were analyzed. The protein–ligand complex mean deviation in atom displacement during a specific time interval from the original (native structure) frame is known as the root mean square deviation (RMSD). Determination of the stability of the protein and ligand structure during simulation is a crucial quantitative assessment. It implies that the protein–ligand combination is more stable if there is less variation in the RMSD during the MD simulation and vice versa (Liu et al., 2018). Protein C-atoms were found to have mean RMSD values of 1.332 Å, 1.461 Å, 1.335 Å and 1.511 Å when bound with ligands referenced in Fig. 9A, ZH-1, ZH-2, and ZH-3; in contrast, the mean RMSD of 1HIV protein C-atoms was determined to be 1.41 Å (Table 9). The ZH-3-1HIV complex exhibited slight variations at 80 ns with an RMSD of 1.4 Å, according to the protein's RMSD analysis (Fig. 9B), before converging to a stable system at the conclusion of the simulation period. Fig. 9C demonstrated that all of the systems (reference, ZH-1, ZH-2, and ZH-3) maintained the molecular interaction profile during the 100 ns MD simulation run. More significantly, they did not exhibit any significant fluctuation, which is unquestionably acceptable to assess the consistency of the protein–ligand complex interaction. Out of all the compounds that were found, reference had the highest RMSD deviations when ligand–protein RMSD values were analyzed. The RMSD in ZH-1 and ZH-2 remained constant during the MD simulation. In ZH-3, equilibration caused initial fluctuations that eventually became constant with only slight variations. The RMSD data suggests that the ligand under research have a strong stabilizing effect on the 1HIV complex protein when compared to the reference (Fig. 9B). The stability of protein–ligand complex analyzed using the RMSF technique. It is a measurement of the flexibility of a protein residue in terms of fluctuations in protein atoms across a 100 ns MD simulation. The average RMSF values for 1HIV bound with reference, ZH-1, ZH-2 and ZH-3 molecules were measured to be 0.953 Å, 0.866 Å, 0.862 Å, and 0.930 Å respectively. However, it was shown that the RMSF values of Apo protein was 0.903 Å (Table 1). It is evident from the low RMSF values that no significant alterations in amino acid residues emerged during the experiment. The ligand made interaction with Asp25, Gly27, Ala28, Gly49, Ile50, and Gly48 in the ZH-1HIV protein complex. There were variations seen in the residues that make up the binding pocket of up to 2.6Å. The reference, ZH-1, ZH-2 and ZH-3 have the least amount of variations. However, reference has a comparable RMSF trend that is consistent with its RMSD pattern. The slight variations in RMSF made it evident that the ZH-1HIV amino acids did, in fact vary considerably at Gly16 and Gly52 throughout the simulation (Fig. 2C). Another criterion for figuring out the root mean square distance of a collection of atoms from their common centre of mass is the evaluation of radius of gyration (Rg), which may provide an indication of the approximate compactness of system. In Fig. 9D, the Rg data for each frame were plotted against time to show how they varied over time. Table 1 displays the maximum, lowest and median Rg values. Fig. 2D shows that for all of the complexes, there was no significant change in Rg. The small differences in Rg values demonstrated that all systems adopt compact protein structures equally, which may be attributable in part to significant contributions from the protein structure α-helices. It was evident from the RMSD, RMSF, and Rg parameters derived from the MD simulation trajectories that 1HIV inhibitors were retained within the protein cavity and that protein–ligand complexes were stable in a dynamic state (Table 9). To gain a better understanding of the inhibitors' binding affinities with 1HIV protein, Throughout the 100 ns simulation, each ligand–protein complex interaction was examined. The protein–ligand interactions histogram for the full 100 ns MD experiment was shown in Fig. 10. The protein–ligand contacts showed that the identified inhibitors primarily interact with HIV through hydrophobic and hydrogen bond interactions. As a result, the binding free energies computation utilizing MMGB/SA calculations further confirmed the observed results.

Close

Fig. 9

(A) Structural changes at C-α atoms of protein (RMSD) over entire run of 100 ns trajectory on binding with identified ligands and reference molecule, (B) RMSD values with protein backbone (pdb id: 1HIV), (C) Fluctuations of residues (RMSF) on binding with ligand, and (D) The compactness of protein represented by radius of gyration.

Export to PPT

Close

Fig. 10

Protein ligands Interactions and protein ligand contacts histogram of (A) reference molecules, (B) ZH-1, (C) ZH-2 and (D) ZH-3 respectively at the active site of HIV over the entire run of 100 ns (PDB id:). The green, violet and blue color represent the hydrophobic interactions, H bond and water bridges respectively.

Export to PPT

3.7 PCA analysis

PCA technique is widely used powerful dimensionality reduction method that can be used to extract useful information from MD simulation data on the main collective movements and conformational changes of biomolecular systems. Large-scale simulations of Molecular Dynamics provide a huge amount of data. To study a medium-sized protein accurately, the positions of at least 10,000 atoms must be calculated every 10⁻¹⁵ s. There are five steps involved in PCA including data centring, covariance matrix computing, eigenvector and eigenvalue computation, eigenvector sorting by eigenvalue sorting and data projection onto the new feature space. The Fig. 11 displayed the PCA analysis of trajectories including reference and ZH1, ZH2, ZH4. The PCA analysis of 1HIV-reference complex showed that PC1, PC2 and PC3 contributed for the 79.1 %, 7.81 %, and 4.01 % respectively (Fig. 11A). Similarly, 1HIV-ZH1 displayed the 78.62 %, 11.92 %, and 3.37 % values for PC1, PC2, and PC3, respectively (Fig. 11B). For 1HIV-ZH2, PC1, PC2, and PC3 values were 74.67 %, 15 %, and 3.51 %, respectively (Fig. 11C). Together, these three top PCs accounted for most of the variation in the original distribution, yielding a succinct yet useful picture of the atomic fluctuations among the studied conformations. As a result, we reduced the complexity of the paths of the simulated systems onto this two-dimensional landscape by concentrating our research on the conformational changes projected onto the subspace defined by PC1 and PC2. The projection allowed the Cα atoms three-dimensional fluctuations to be captured in a lower-dimensional form, which made it easier to understand the conformational shifts. The continuous colour scale, which went from blue to white to red, showed irregular transitions between different conformational states. A clear difference between the blue and red conformations could be seen by projecting the trajectories along the PC1 and PC2 axes which suggests that the simulated systems underwent a conformational transition. The conformational diversity sampled during the simulations was brought to light by this approach, which also offered insights into possible conformational routes and transitions connected to ligand binding and unbinding events.

Close

Fig. 11

Instantaneous conformations of trajectory obtained by PCA analysis coloured in chronological order from blue to red. The degree of conformational changes in the simulation was described by the spread of blue and red color dots where the color ranges from blue to white to red corresponds to simulation duration. The first, intermediate and the final timestep is represented by blue, white and red respectively. The atomic displacements for PCA-1 are shown on the left side, indicating notable changes. (A) A reference molecules, (B) ZH-1, (C) ZH-2 and (D) ZH-3.

Export to PPT

3.8 Binding free energies calculations (MMGB/SA)

Molecular mechanics/Generalized Born Surface Area (MMGB/SA) is a critical computational technique in drug design and molecular modelling. It is essential to studying and maximizing molecular interactions, particularly in the context of drug discovery. MMGBSA is an essential tool for rational drug design and discovery projects because it provides quantitative insights into molecular interactions, helps in the identification of lead compounds, improves their binding affinities and guides decision-making throughout the drug development pipeline. It contains contributions from molecular mechanics (MM) interactions and solvation effects (Genheden and Ryde, 2015). The calculated dG Bind scores were found to be −64.2064, −100.314, −86.211, and −100.461 kcal/mole, respectively for reference, ZH-1, ZH-2 & ZH-3 molecules. The least binding energy for ZH-3 (−100.461 kcal/mole) may be due to the additional interactions what we have discussed in MD simulation study part. All of the identified HITs exhibit good interactions, according to the additional binding energies such as dG_Bind_Coulomb, dG_Bind_Covalent, dG_Bind_Hbond, dG_Bind_Lipo, and dG_Bind_Packing (stacking interactions) which were also discovered in negative values. All the values were summarized in Table 10. The substances with lower (more negative) total binding free energies are great candidates for further investigation because they are more likely to interact strongly with the target protein.

In-silico ADME calculation: In this study, in-silico ADME prediction of Reference, ZH-1, ZH-2 and ZH-3 was done with the results summarized in Table 11. The parameters mentioned in material and methods part have been evaluated. Additionally, the newly synthesized compounds did not cross the blood–brain barrier and were not predicted to have central nervous system toxicity. For the most active predicted compound ZH-1 was found to have a 0 & 6.5H bond acceptor and donor respectively. The obtained results revealed that molecule has a lower oral absorption rate and metabolism. Moreover, ZH-1 contained the 83.459 polar surface area. Apart from this other physicochemical properties QPlogPo/w 6.3, QPlogBB of −0.761 and a QPPMDCK of, suggesting identified molecule has a minimal toxicity and all drug likeness properties.