3.1 Screening of drug molecules based on pharmacokinetics
The SuperNatural library contains 449,058 natural compounds along with their structural and physicochemical information[34]. Based on Lipinski's rule of five, 21104 drug-like molecules were screened. These molecules, along with 800 NatProd molecules, were subjected to further screening by predicting their ADMET properties. Most of the compounds qualified for the Golden Triangle rule, Lipinski’s rule of five, and Pfizer’s rule. Molecules with predicted carcinogenicity, rat oral acute toxicity, AMES toxicity, and mutagenicity potentials were removed, and finally, 1,104 molecules were selected for structure-based virtual screening (Table S1).
3.2 Structure prediction by homology modeling
In the model quality assessment, the predicted structures of DNA polymerase β from the amino acid sequences showed > 96% favorable regions in the Ramachandran plots, with the QMEAN score < 0.90 (> 0.6), as determined by the MolProbility tool of the SWISS MODEL.
3.3 Pocket Identification and Docking-Based Virtual Screening
Grid-based HECOMi finder (Ghecom)is a program for finding multi-scale pockets on protein surfaces using mathematical morphology[34, 54, 55]. A total of fifteen pockets were identified in the structure of the DNA pol β protein by the Ghecom algorithm. The largest pocket in the DNA pol β protein had a volume of 3654.65 Å3. The second largest pocket had a volume of 246.78 Å3. The smallest pocket within the structure of the pol β protein had a volume of 29.18 Å3. In this study, the second largest pocket was chosen which lies into the catalytic domain of the protein, were set as (X:3.672; Y:12.247; Z:6.153).
After the primary screening, approximately 1104 drug molecules were selected. The selected molecules were screened using docking-based virtual screening using the EsayVS server. Compounds with an affinity score of -6.5 with h-bond 7 protein-ligand interactions were selected. Totals of 135 molecules were finally selected against pol β protein (Table S2). These 135 molecules were further screened using the DockThor server (Table S3).
Additionally, the selected molecules were further screened based on their EC using Flare v5.0.0, and compounds with EC scores > 0.25 were considered for further study.
Table 2
Molecular Docking Parameters and Binding Characteristics.
ID
|
Easy Vs
|
ThorDock
|
flare
|
Affinity
|
H bonds
|
Internal energy
|
Total energy
|
Van der waal energy
|
Electrostatic energy
|
EC score
|
Amino acid involved in h bond
|
No of π-π bond
|
Aa in π-π bonds (bond distance)
|
SN00305418
|
-7.6
|
8
|
-39.107
|
-1.302
|
-9.785
|
-29.322
|
0.349
|
SER (A204), LYS (A206), LYS (A234), THR (A233), LEU (A259)
|
2
|
PHE(A200)(2.8), LYS(A234)(3.0)
|
SN00158342
|
-7.6
|
8
|
-35.5
|
-10.752
|
-4.305
|
-31.195
|
0.333
|
SER (A202), LYS (A206), GLU (A295)
|
2
|
LYS(A206)(2.6), TYR(A296)(3.1)
|
SN00004251
|
-7.1
|
12
|
-42.452
|
8.728
|
-16.872
|
-25.58
|
0.335
|
LYS (A206), GLU (A232), LEU (A259), SER (A204), THR (A201)
|
2
|
THY(A233) (3.7), PHE(A200)(3.0)
|
SN00341636
|
-7.4
|
9
|
-43.26
|
24.963
|
-20.089
|
-23.171
|
0.294
|
ASP (160), GLU (A329), VAL (A177), ARG (A152), ARG (A182)
|
0
|
|
3.4 selection of drug molecules based on Neural Network
Based on 6 parameters, affinity score, internal energy, van der Waal force, EC scores, electrostatic interaction, and no of h bonds, finally, seven molecules, SN00261400, SN00305418, SN00006989, SN00158342, SN00305418, SN00004251, and SN00341636 were selected from the pool of 135 compounds (Table 3). As SN00006989 and SN00158342 were found to be structurally similar to the SN00261400, and SN00305418 was found to be structurally similar to the SN00305418, only four molecules, SN00305418, SN00158342, SN00004251, and SN00341636 were considered for the md simulation study.
Table 3
Selected drug molecules based on Neural Network
No.
|
1
|
2
|
3
|
4
|
Supernatural Id
|
SN00305418
|
SN00158342
|
SN00004251
|
SN00341636
|
3.5 Analysis of Protein-Ligand Interactions
The intermolecular interactions between the selected molecules with the binding sites of the Pol β are analyzed from the best docked pose for each molecule.
The SN00158342 molecule is primarily stabilized by four hydrogen bonds to Ser202, Glu295, and Lys206, with bond distances ranging from 1.9 to 2.3 Å. Additionally, it forms two van der Waals bonds with Lys206 and Tyr296, having bond distances of 2.6 Å and 3.1 Å, respectively.
Similarly, SN00305418 is stabilized by six hydrogen bonds and two van der Waals bonds. The molecule forms hydrogen bonds with Ser204, Lys206, Thr233, Lys234, and Leu259, with bond distances ranging from 1.8 Å to 2.5 Å. It also creates two van der Waals bonds with Phe200 and Lys234, having bond distances of 2.8 Å and 3.0 Å, respectively.
The SN00004251 forms five hydrogen bonds with Thr201, Ser204, Glu232, and Leu259, with bond lengths ranging from 1.9 to 2.2 Å. Additionally, two van der Waals bonds are observed between the molecule and Phe200, and Thr233, with bond distances of 3.0 Å and 3.7 Å, respectively.
Unlike the previous three molecules, SN00341636 does not form van der Waals bonds with the target protein. Instead, the ligand is stabilized by forming five hydrogen bonds with Arg152, Val177, Arg182, Glu329, and Pro330, with bond distances ranging from 1.9 Å to 2.4Å.
3.6 MD Simulations of Protein-Ligand Complexes
MD simulation is a computational approach to predict and analyze the stabilities of protein-ligand complexes, and to study the atomic movements with respect to a macromolecule. The stabilities and behaviors of the protein-ligand complexes were analyzed in a dynamic environment based on the following parameters; (i) root-mean-square deviation (RMSD), (ii) root-mean-square fluctuation (RMSF), (iii) the radius of gyration (Rg), and (iv) molecular mechanics/generalized Born surface area (MM/GBSA) energy[52, 55].
3.6.1 RMSD Values of the Cα Backbone of Target Proteins
The RMSD values of the four protein-ligand complexes were visualized by plotting RMSD values against time, which helps to understand the structural stability and integrity of the complexes. The trajectory analysis revealed that all the four compounds complex with pol β remain stable throughout the 100 ns simulation. The average RMSD values of the 4 molecules selected against pol β ranged from 2.449 to 4.781 Å. Compounds SN00004251, SN00341636 possess high stabilities during the simulation, with an RMSD fluctuation of 0.378 Å and 0.444 Å respectively. (TableS5).
The SN00305418, SN00158342 molecules found to be less stable during the simulation, with an RMSD fluctuation of 0.837 Å and 0.602 Å respectively. The overall RMSD fluctuation of protein is more in comparison with other four protein ligand complexes.
3.6.2 RgValues of the Cα Backbone of Target Proteins
The compactness of the protein-ligand complexes during the simulation was determined by measuring the values of Rg. The average values of Rg for the 4 compounds complexed with pol β ranged from 22.177 to 22.636 Å (Table S6). Rg of pol β fluctuates more in comparison with other four protein ligand complexes. Suddenly the Rg value of pol β increase after 80 ns of simulation.
3.6.3 RMSF Values of the Cα Backbone of Target Proteins
The average atomic mobility of the protein backbone during the MD simulations was measured using the values of RMSF. The average RMSF values of the 4 molecules complexed with pol β ranged from 1.703 to 2.020 Å. Further analysis revealed that residues 203, 246 and 306 of pol β underwent fluctuations(Table S7).
However, only two amino acid residues 246, and 306 underwent fluctuations when complexed with SN00305418, SN00158342, SN00004251 and SN00341636. Residues Arg152, Asp160, Val 177, Thr201 Ser204, Lys206, Lys234, Thr233, Leu259, Glu295 and Glu329 which were primarily involved in the formations of ligand–protein hydrogen bonds with four different chemical compounds, underwent minimal fluctuations. Similarly, the amino acid residues that mediated the formations of stacking interactions remained stable during the simulation.
3.7 Determination of Binding Free Energies of Protein-Ligand Complexes
The binding free energy represents the sum total of all the interaction energies, including the van der Waals energy, polar solvation energy, electrostatic energy, and solvent accessible surface area SASA energy. The binding free energies of all the complexes were estimated using the MM/GBSA approach. The binding free energies of the four compounds complexed with pol β ranged from − 15.9181 to -29.7465 kcal/mol, of which SN00004251 (-29.7465 ± 6.7833) had the lowest free energy of binding.