2.1. Sequence Alignments and Molecular Modeling
Among the BLASTp results, the structure was selected as templates: homoserine dehydrogenase from Thiobacillus denitrificans
(PDB code 3MJT). The sequence identities between HSD and templates 3MJT was 40%. It well known that above 50% sequence identity, models tend to be reliable, with only minor errors in side chain packing and rotameric state. In the 30%–50% identity range, errors can be more severe and are often located in loops. Below 30% identity, serious errors occur, sometimes resulting in the basic fold being mis-predicted [18
]. Thus, 40% identity is sufficient homology to construct a believable model [19
]. The sequence alignment performed using the MolsoftICM for homology modeling is shown in Figure 1a
. Although the sequence identities between HSD and templates: Putative Homoserine Dehydrogenase (NP_069768.1) from Archaeoglobus Fulgidus
(PDB ID 3DO5) (41%) is higher than that of HSD and 3MJT’s, 3MJT was chosen as template. The reasons are as follows: firstly, 3MJT contains 496 residues, and 3DO5 contains 327 residues. As a template, the length of 3MJT is more appropriate than that of 3DO5’s. Secondly, phylogenetic analysis (seen from Figure 1b
) showed that 3MJT and HSD are the same subfamily, and thus their spatial structure should be more similar. The most significant step in homology modeling process is to obtain the correct sequence alignment of the target sequence with the homologues, and it reveals that the residues involved in binding of substrate in templates (Lys217 (Proton donor)), Arg117 (NADH binding) and Glu196 (substrate binding site) were conserved (the corresponding residue: Lys207, Arg107 and Glu192) in HSD.
The coordinates of the crystal structures of homoserine dehydrogenase from Thiobacillus denitrificans
(PDB code 3MJT) was used as templates to build the structure of HSD. The 3D model of the HSD was built by Swiss model [16
]. Further, refinement was performed in order to obtain the best conformation of the developed model of HSD. Analysis of 20 ns dynamics shows that the HSD structure is stable and indicated that the system is stable.
The superposition of the average structure of the HSD with the initial model Figure 2
. Figure 2
does not show major structure conformational changes in comparison to the initial model, which is consistent with the relatively low RMSD values. We selected the average structure of the HSD through the further study.
2.2. Validation of Homology Model
The first validation was carried out using Ramachandran plot calculations computed with Molprobity program by checking the detailed residue-by-residue stereo-chemical quality of a protein structure [20
]. The results are shown in Figure 3
. Altogether, 95.0% of all residues were in favored regions, and 98.4% of all residues were in allowed regions. In comparison with the homology model, the template, 3MJT, had a similar Ramachandran plot 98.35% in the allowed regions.
Seen from Figure 3
, it can be concluded that Glu263, Leu62, Tyr340, and Gln264 are in the disallowed regions. So we fixed them by hand. ERRAT is a so-called “overall quality factor” for nonbonded atomic interactions, and higher scores mean higher quality [23
]. The normally accepted range is >50 for a high quality model [23
]. In the current case, the ERRAT score for HSD model is 87.50, well within the range of a high quality model, in the mean time the ERRAT score for the templates 3MTJ is 96.95. Thus, the above analysis suggests that the backbone conformation and non-bonded interactions of HSD homology model is all reasonable within a normal range. The final evaluation of the built HSD structure was checked by Verify 3D [24
represents the Verify 3D graph of the predicted HSD model. It is to be noted that compatibility scores above zero correspond to acceptable side chain environment. From Figure 4
, we can see that almost all residues are reasonable. In brief, the geometric quality of the backbone conformation, the residue interaction, the residue contact and the energy profile of the structure is all well within the limits established for reliable structures. All evaluations suggest that a reasonable homology model for HSD has been obtained that can be exposed for examination of protein-substrate and protein inhibitor interactions.
2.3. Identification of Substrate-Binding Region and Co-Factor-Binding Region in HSD
HSDs have a common co-factor, NADH. It was reported that the NADH is bound to the Rossmann fold in the conventional mode, that is, the cofactor–enzyme interactions are predominantly mediated through hydrogen bonds between cofactor phosphate moieties and sugar hydroxyl groups with enzyme amide backbone groups [4
]. However, the orientation of NAD+
of the nicotinamide ring is consistent with the experimentally observed stereospecificity of hydride transfer, that is to say that NADH hydride is facing the substrate-binding region.
The cavity volume estimated by CASTp [26
] is dependent on the radius of the probe sphere; a probe radius of 1.4 Å outlines a cavity of 1124.2 Å3
for HSD of substrate binding, while a probe radius of 1.4 Å outlines a cavity of 1367.2 Å3
for 3MJT of substrate binding. From this result, we can conjecture that the active site of HSD is almost the same size of 3MJT.
Compared with the residues in the NADH binding site of the other HSDs’, the residues of HSD participating in the NADH binding site are listed in Table 1
. In order to confirm whether the binding site determination for NADH is correct, the residues in the NADH-binding site of homoserine dehydrogenase (PDB ID 1EBF) [4
] were coordinated with NADH that have been determined by X-ray. Seen from Table 1
, we can conclude that the corresponding residues in HSD are Thr313, Val18, Gly19, Gly45, Gly47, Ile13, Ile46, Gly16, Gly14, Asn17, Ala105, Asn106, and Arg111. Half of these residues are conserved, and the reason may lie in the low sequence identity (27%). Hence the different NAD+
-binding site residue between 1EBF and HSD may affect the catalytic efficiency of the substrate of these two enzymes.
2.4. Docking Study
The substrate, Aspartate-4-semialdehyde (ASA), and NADH are docked to HSD with Autodock 4.2 [27
] and AutoDock vina [29
], respectively. The grid size for Autodock 4.2 is 56 × 56 × 56 Å, and grid size for AutoDock vina is 24 × 24 × 24 Å. The results are listed in Table 2
Seen from Table 2
, the docking score is −5.24 Kcal·mol−1
for Autodock 4.2, while for AutoDock vina the docking score is 4.50 Kcal·mol−1
. The RMSD between the original and the docked of the 3D structure of ASA are also shown in Table 2
(27.10 Å for AutoDock vina, 25.20 Å for Autodock 4.2). These results showed the binding mode generated by Autodock 4.2 is more reasonable than that of AutoDock vina, and thus chosen for further study.
shows the substrate, ASA, and NADH docking in the HSD. We can see that ASA and NADH are in the groove. ASA is located near NADH, and so it is useful to transfer action.
Hydrogen bonds may be important in substrate binding. There are four hydrogen bonds between ASA and HSD (seen from Figure 6
and Table 3
). Thr163 make two hydrogen bonds with the ASA (1.97 and 1.70 Å). Asp198 forms a strong hydrogen bon with the NH group with ASA. There is a weak hydrogen bond (2.26 Å) between the NH group of Glu192 and O8 atom of ASA. From Figure 6
, it can be seen that Ala191 and Ala308 have electronic contact with ASA, whereas Tyr190, Glu192, Asp198, Gly162, Thr163, and Asn161 have strong van der Waals (VDW) contact with ASA. So the binding pocked of ASA to HSD may contain Ala191 and Ala308, Tyr190, Glu192, Asp198, Gly162, Thr163, and Asn161. In particular, Thr163, Asp198, and Glu192 may be important for ASA binding for they make hydrogen bonds with HSD. Glu192 is substrate binding residue, and this result is completely consistent with experimental data [15
2.5. High Throughput Virtual Screening Procedure and Docking the Inhibitor to the Protein
It was reported that the receiving operating characteristic (ROC) evaluation is described as the ratio of the true positive rate to the false positive rate when a given proportion of known decoys have been observed [31
]. Moreover, the receiving operating characteristic (ROC) curve, a graphical plot of the sensitivity (true positive rate, sensitivity) VS specificity (false positive rate), was calculated to avoid the sensitivity for small changes in ranking, where T score represents the number of correctly identified actives, and F score (false positives) represents the number of decoys incorrectly predicted as actives. The area under the ROC curve (AUC plot) gives the probability of ranking a randomly selected active higher than a randomly chosen decoy [33
]. It ranges from 0 to 1, where 1 indicates a perfect ranking, where in all actives ranked above the decoys. A traditional academic point system for classifying the accuracy of a virtual screening test is known as follows: AUC < 0.5 is fail; 0.5 ≤ AUC < 0.70 is poor; 0.7 ≤ AUC < 0.8 is fair; 0.8 ≤ AUC < 0.9 is good; and 0.9 ≤ AUC ≤ 1 is excellent [34
]. Essentially independent of the actual number of positive and negative instances, the AUC of an ROC plot gives an objective measure of query performance. In this study, 32 inhibitors [14
] were used to generate the decoys using DUD-E on line [31
]. Autodock 4.2 [27
], Autodock vina [29
] and Dock 3.6 [37
] are used for docking. The ROC curve is shown in Figure 7
. From Figure 7
, AUC plot with Autodock 4.2 is 0.64, which is larger than that of AutoDock vina and Dock 3.6. And so Autodock 4.2 is used to further virtual screening.
Virtual screening of compound libraries has become a standard technology in modern drug discovery pipelines. The “2008/5” version of a Natural Products Database (NPD) contains almost 90,000 commercially available compounds. The target used in our study was the 3D structure of HSD mentioned above. In this simulation, it was also been screened that the potent inhibitor of HSD, 4-(4-HYDROXY-3-ISOPROPYLPHENYLTHIO)-2-ISOPROPYLPHENOL (178), which is found to be competitive with ASA (IC50
5.1 μm) [15
], as the leader drugs searching in Zinc data for 50% similarity.
After the screening, 164 compounds have been found. AutoDock 4.2 is used for virtual screening. AutoDock 4.2 uses a semi-empirical free energy force field to evaluate conformations during docking simulations. The force field was parameterized using a large number of protein-inhibitor complexes for which both structure and inhibition constants, or Ki
are known. Table 4
listed the discovered inhibitors from the docking screen against the known crystal structure [14
]. Two inhibitors showed good inhibition with the substrate ASA [14
]. This result was consistent with our docking results (the free energy of binding and the calculated Ki
) (Table 4
shows the binding pose of the inhibitor 178 in the HSD. In particular, Lys107 has a cation–π interaction with the inhibitor 178. The flat face of an aromatic ring has a partial negative charge owing to the pi electrons. Cations such as the sidechains of Lys or Arg, cationic ligands, or metal cations often align themselves centered over the faces of aromatic rings. It was reported that cation–π interactions should be considered alongside the more conventional hydrogen bonds, salt bridges, and hydrophobic effects in any analysis of protein structure, and they can also contribute significantly to intermolecular contacts and interactions with ligands [38
]. As discussed before, the cation–π interaction makes the inhibitor-enzyme stable and cannot be removed.
As shown in the Figure 9
, in the HSD-178 complex, Ala105, Lys107, Glu192, Thr163, Asp198, and Tyr191 have electronic contact with inhibitor 178, and Asn161, Gly162, K207, and Asp203 have strong VDW contacts with inhibitor 178. In particularly, Lys107 forms a cation–π interaction with the inhibitor 178, and this result indicates Lys107 is an important residue for inhibitor 178. These results can serve as a guide to the selection of candidate sites for further experimental studies of site directed mutagenesis.
lists four compounds. The binding energies and calculated Ki
between four compounds and HSD are all lower than that of 178-HSD’s (−6.07 Kcal·mol−1
, 17.56 μM). The similarity of the new compounds versus
178 was assessed by calculating the Tanimoto coefficient (Tc
) with Discovery studio 3.5 client to the 164 HSD inhibitors annotations. Tc
value ranges from 0 to 1, where 0 represents no detection of the same bits; however, 1 does not mean that the two molecules are totally identical. The atom pair similarities (SimAB
) will be determined by the number of atom pair types shared by the two molecules, where 0 indicates no similarity and 1 indicates identity [32
]. The four top-scoring docking hits (Tc
< 0.6) is selected. The four top-scoring docking hits all seemed to cation–π ion pair with the key recognition residue Lys107, and Lys207 (from Figure 9
). These ligands therefore seemed to be new chemotypes for HSD.