Discovery of Novel DUSP16 Phosphatase Inhibitors through Virtual Screening with Homology Modeled Protein Structure

Abstract

Recently, dual-specificity phosphatase 16 (DUSP16) emerged as a promising therapeutic target protein for the development of anti-atherosclerosis and anticancer medicines. The present study was undertaken to identify the novel inhibitors of DUSP16 based on the structure-based virtual screening. We have been able to find seven novel inhibitors of DUSP16 through the drug design protocol involving homology modeling of the target protein, docking simulations between DUSP16 and its putative inhibitors with the modified scoring function, and in vitro enzyme assay. These inhibitors revealed good potency, with IC₅₀ values ranging from 1 to 22 µM, and they were also screened computationally for having desirable physicochemical properties as drug candidates. Therefore, they deserve consideration for further development by structure-activity relationship studies to optimize the inhibitory activity against DUSP16. Structural features relevant to the stabilization of the newly identified inhibitors in the active site of DUSP16 are addressed in detail.

Keywords

cheminformatics computational chemistry enzyme assays medicinal chemistry

Introduction

The mitogen-activated protein kinase (MAPK) pathway controls the cellular functions associated with the growth, differentiation, and death of cells.¹ In this signaling cascade, MAPKs are activated by the phosphorylation in the active loop residues. MAPK phosphatases (MKPs) dephosphorylate MAPKs and thereby turn off their activities.² Of various MKPs, dual-specificity phosphatase 16 (DUSP16; also known as MKP-7) dephosphorylates JUN kinases (JNKs) specifically and thereby regulates T-cell differentiation,³ which can be responsible for the pathogenesis of Burkitt’s lymphoma⁴ and cerebral ischemia.⁵ Recently, DUSP16 was also found to induce the expression of VCAM-1 through the suppression of JNK signaling.⁶ Because VCAM-1 plays a pivotal role in the development of atherosclerosis and tumor progression,⁷ DUSP16 can serve as a potential target for the development of anti-atherosclerosis and anticancer medicines. The merit of DUSP16 as a target for drug discovery was further supported by the experimental finding that the knockout of DUSP16 gene by small interfering RNA inhibited the expression of VCAM-1.⁶

Although the three-dimensional (3D) structure of DUSP16 is not yet known, X-ray crystal structures of its highly homologous proteins such as DUSP9 were reported in complex with a substrate analogue.⁸ To accommodate the small phosphorylated substrate group, DUSP9 has a planar and shallow active site in which the negatively charged substrate group can be stabilized by the hydrogen bonds with the backbone amide groups and the side-chain guanidinium ion of the Arg residue located near the catalytic cysteine residue (Cys290). Structural information about the active site and the interactions with a small-molecule ligand is very useful for designing new potent inhibitors that can develop into drug candidates. Indeed, many potent and selective small-molecule phosphatase inhibitors have been identified based on such a structure-based molecular design approach.^9,10 Nonetheless, the discovery of DUSP16 inhibitors has lagged behind the biological and pharmaceutical studies. No small-molecule DUSP16 inhibitor has been reported so far in the literature. This may be attributed in a large to the lack of 3D structure of DUSP16, which makes it difficult to apply the rational design approaches for finding the novel inhibitors.

In this study, we try to identify the potent DUSP16 inhibitors through the structure-based drug design protocol, which involves homology modeling of the 3D structure of DUSP16, virtual screening of a chemical library with docking simulations, and an in vitro enzyme assay. Virtual screening using molecular docking has not always been useful because of the inaccuracy in the scoring function that estimates the binding affinity between the target protein and a putative inhibitor.¹¹ Therefore, prior to carrying out the virtual screening, we improve the existing scoring function by the implementation of an accurate solvation model to reflect the effect of desolvation cost for a ligand on the binding affinity to the target protein. This modification seems to have an effect of increasing the hit rate in the enzyme assay because the neglect of the solvation term often leads to the overestimation of the binding affinity of a ligand with many polar atoms.¹² We find in this study that the docking simulations with the improved scoring function can be a useful computational tool for elucidating the activities of the identified inhibitors, as well as for enriching the chemical library with molecules that are likely to have desired biological activities.

Materials and Methods

Because of the lack of structural information on DUSP16, we performed homology modeling using the X-ray crystal structure of DUSP9 (PDB entry: 2HXP) [8] as the template to build the 3D structure of DUSP16 necessary for the structure-based virtual screening. The sequence alignment between the phosphatase domains of DUSP9 and DUSP16 was derived with version 2.1 of the ClustalW program,¹³ using the BLOSUM matrices for scoring the alignments. Opening and extension gap penalties were changed systematically, and the derived alignments were inspected for the violation of structural integrity in the structurally conserved regions. More specifically, we adopted only the alignments in which the sequences of active-site residues of DUSP9 (CLAGVSR) and DUSP16 (CLAGISR) were matched because all DUSPs share the catalytic core in a form of CXXGXXR. Using the best-scored sequence alignment, the 3D structure of the catalytic domain of DUSP16 was constructed using the MODELLER program (version 9.12).¹⁴ In this homology modeling, we applied the optimization method involving conjugate gradients and molecular dynamic simulations to minimize the violations of the spatial restraints. With respect to the structure of the gap regions, the coordinates were built from a randomized distorted structure that resides approximately between the two anchoring regions as implemented in the MODELLER program. To refine the calculated structures, loop modeling was also carried out based on the enumeration algorithm.¹⁵ WPD (residues 209–216) and Q (residues 279–285) loops were included in this structural refinement. As the final structural model for DUSP16, we selected the one with the lowest value of MODELLER objective function. Finally, we calculated the conformational energy of the predicted structure of DUSP16 with the ProSa 2003 program¹⁶ for the purpose of evaluation.

The atomic coordinates obtained from the homology modeling were used as the receptor model in the virtual screening of DUSP16 inhibitors. Special attention was paid to the determination of the protonation states of the ionizable residues of the homology-modeled structure. The side chains of Asp and Glu residues were assumed to be neutral if either of the carboxylate oxygens was directed toward a hydrogen-bond accepting group such as the backbone aminocarbonyl oxygen within the distance of 3.5 Å, which is a generally accepted distance limit for the hydrogen bond of moderate strength.¹⁷ Similarly, the side chains of lysine residues were assumed to be protonated unless the NZ atom stayed in proximity to a hydrogen-bond donating group. The same procedure was also applied for determining the protonation states of His residues. As a consequence of the inspection for protonation states, Asp199, Lys234, and Asp264 were found to reside in the neutral form in the final structural model of DUSP16. After the assignment of the protonation states, we carried out 200 cycles of energy minimization for DUSP16 with the AMBER program of version 12 to remove the bad steric contacts. In this step, it should be kept in mind that the energy minimization should be terminated once all bad van der Waals contacts were removed because the full minimization in the gas phase can make the protein structure converge to a physically unacceptable one because of the overestimation of intramolecular electrostatic interactions.

The docking library for DUSP16 comprising about 260,000 compounds was constructed from the chemical database distributed by Interbioscreen (http://www.ibscreen.com) containing approximately 477,000 synthetic and natural compounds. Prior to performing the docking simulations for virtual screening, they were filtrated according to Lipinski’s “Rule of Five” to select only the compounds with the physicochemical properties of potential drug candidates¹⁸ and without reactive functional group(s). To remove the structural redundancies in the chemical library, the structurally similar compounds with Tanimoto coefficient larger than 0.8 were clustered into a single representative molecule. As a consequence of the filtering process, a docking library containing ~260,000 compounds was constructed. These compounds were then processed with the CORINA program to obtain their 3D atomic coordinates, followed by the assignment of atomic charges with the Gasteiger-Marsilli method.¹⁹ We used version 3.5 of the AutoDock program²⁰ in the virtual screening of DUSP16 inhibitors. AMBER force field parameters were used for calculating the van der Waals interactions and the internal energy of a ligand as implemented in the original AutoDock program. Docking simulations with AutoDock were then carried out in the active site of DUSP16 to score and rank the compounds in the docking library according to their calculated binding affinities.

In the docking simulations for virtual screening, we used the modified AutoDock scoring function to which the solvation free energy term for organic molecules was added to estimate the desolvation cost for protein-ligand association. This modified scoring function can be written in the following form.

\begin{array}{l} Δ G_{b i n d}^{a q} = W_{v d W} \sum_{i = 1} \sum_{j = 1} (\frac{A_{i j}}{r_{i j}^{12}} - \frac{B_{i j}}{r_{i j}^{6}}) + W_{h b o n d} \sum_{i = 1} \sum_{j = 1} E (t) (\frac{C_{i j}}{r_{i j}^{12}} - \frac{D_{i j}}{r_{i j}^{10}}) \\ + W_{e l e c} \sum_{i = 1} \sum_{j = 1} \frac{q_{i} q_{j}}{ε (r_{i j}) r_{i j}} + W_{t o r} N_{t o r} + W_{s o l} \sum_{i = 1} S_{i} (O_{i}^{\max} - \sum_{j \neq i} V_{j} e^{- \frac{r_{i j}^{2}}{2 σ^{2}}}) \end{array}

Here, W_vdW, W_hbond, W_elec, W_tor, and W_sol are the weighting factors of van der Waals, hydrogen bond, electrostatic interactions, torsional term, and solvation free energy of a putative inhibitor, respectively. r_ij represents the interatomic distance, and A_ij, B_ij, C_ij, and D_ij are associated with the depths of the potential energy well and the equilibrium separations between the two atoms. The hydrogen bond term has an additional weighting factor, E(t), to represent the angle-dependent directionality. With respect to the distance-dependent dielectric constant (ϵ(r_ij)), a sigmoidal function proposed by Mehler et al.²¹ was used for the calculation of the interatomic electrostatic interactions between DUSP16 and the putative inhibitors. In the entropic term, N_tor denotes the number of rotatable bonds in the ligand. Although the original AutoDock scoring function included a simple solvation energy term, we replaced it with the improved one shown in the last term of eq 1. In this new solvation free energy term, S_i and V_i are the solvation parameter and the fragmental volume of atom i,²² respectively, whereas O_i^max stands for the maximum atomic occupancy that means the volume of atom i exposed to bulk solvent in the absence of the other solute atoms. The negative and positive values of S_i parameter indicate the stabilization and destabilization of the solute atom i, respectively, due to the combined effects of intermolecular interactions with water molecules and intramolecular interactions with the rest of the solute atoms. In the calculation of the molecular solvation free energy term in eq 1, we used the atomic parameters developed by Choi et al. because they proved to be successful in predicting the solvation free energies of a variety of organic molecules.²³ This modification of the scoring function seems to increase the accuracy in virtual screening of DUSP16 inhibitors because the underestimation of ligand solvation often leads to the overestimation of the binding affinity of an inhibitor with many polar atoms. Indeed, the superiority of this modified scoring function to the previous one was well appreciated in recent studies for virtual screening of kinase and phosphatase inhibitors.^24,25

The catalytic domain of DUSP16 was overexpressed in Escherichia coli and purified. For this purpose, E. coli cells containing the DUSP16 expression vector were induced with 0.2 mM IPTG at 18 °C overnight. Cell pellets were lysed in a buffer containing 50 mM Tris-HCl (pH 7.5), 500 mM NaCl, 1% phenylmethylsulfonyl fluoride, 4 mM 2-mercaptoethanol, and 5% glycerol. The His-tagged DUSP16 was purified by Ni-NTA agarose column and dialyzed against 20 mM Tris-HCl (pH 8.0), 200 mM NaCl, 2 mM DTT, and 5% glycerol.

A total of 148 compounds selected from the virtual screening were evaluated for in vitro inhibitory activity against the purified DUSP16. These phosphatase enzyme assays were carried out by using 6,8-difluoro-4-methyl-umbelliferyl phosphate (DiFMUP) as a fluorogenic substrate, which has a K_m value of approximately 20 µM for protein tyrosine phosphatases (PTPs).²⁶ To determine the conditions for the phosphatase assay, we performed the preliminary assay experiments with varying concentrations of DUSP16 and DiFMUP. Then, we selected a combination of the concentrations (3 µM DUSP16 and 10 µM DiFMUP) that generated an ideal progression curve for fluorescence measurements. Under these assay conditions, the fluorescence intensity showed a linear increase in the early time region and reached the maximum after 10 to 20 min. In case of the reaction buffer, we used the conditions under which most DUSPs were found to exhibit an optimal activity.²⁷ Actually, the usefulness of the assay conditions described above was confirmed by elucidating the potency of various known phosphatase inhibitors such as NSC 95397 and Cpd 5.²⁸

The purified DUSP16, DiFMUP, and a candidate inhibitor were incubated in the reaction mixture containing 20 mM Tris-HCl (pH 8.0), 0.01% Triton X-100, and 5 mM DTT for 20 min. The resulting fluorescence was measured by using the PerkinElmer 2030 instrument with excitation and emission wavelengths of 355 and 460 nm, respectively. Initial velocities of reactions were estimated for various inhibitor concentrations (0–50 µM). The inhibitory activities of the putative DUSP16 inhibitors were measured in duplicate at the concentrations of 0.0, 0.625, 1.25, 2.5, 5, 10, 20, and 50 µM to obtain the dose-response curve fits. The IC₅₀ value of each inhibitor was then determined from direct regression analysis using the four-parameter sigmoidal curve as implemented in the SigmaPlot program.

Results and Discussion

The peptide sequence of human DUSP16 comprising 665 amino acids was retrieved from UniProtKB protein knowledgebase (http://www.uniprot.org, accession number: Q9BY84). Although the 3D structure of DUSP16 has not been known, we found from a database search at NCBI BLAST (http://blast.ncbi.nlm.nih.gov) that DUSP9 could serve as a structural homolog for DUSP16. As shown in the supplementary data (Suppl. Fig. S1), the identity and similarity between the amino acid sequences of the phosphatase domains of DUSP9 and DUSP16 amount to 42.5% and 71.9%, respectively. Such a high sequence identity between the template and the target indicates that the structure of DUSP16 can be predicted with accuracy through homology modeling.²⁹ To obtain the atomic coordinates of DUSP16 required for virtual screening with docking simulations, therefore, we performed the homology modeling of the phosphatase domain of DUSP16 using the X-ray crystal structure of DUSP9 as the template.

The structure of the phosphatase domain of DUSP16 generated from the homology modeling was then evaluated with the ProSa 2003 program¹⁶ by investigating whether the interactions of each amino acid residue with the rest of the protein structure could be maintained favorably. This program has been widely used for estimating the reasonableness of a predicted protein structure because the calculated knowledge-based mean fields could serve as a yardstick to judge the quality of protein folds. Figure 1 compares the ProSa energy profiles of the homology-modeled structure of DUSP16 and the X-ray crystal structure of DUSP9. It is seen that the ProSa energy remains negative for all amino acid residues of the template and the target, indicating that both protein structures would be physically acceptable. We also note that the homology-modeled DUSP16 structure appears to be relatively unstable in the N-terminal region when compared with the X-ray structure of DUSP9, and vice versa in the C-terminal region. The overall similarity of the conformational stability of the target to that of the template indicates that the 3D structure of DUSP16 can be predicted with accuracy from homology modeling.

Figure 1.

Comparative view of the ProSa energy profiles for the X-ray crystal structure of DUSP9 and the homology-modeled DUSP16 structure. For convenience, the amino acids are renumbered from 1 instead of retaining the original numbers.

Although it is known that the accuracy of virtual screening can be enhanced significantly by taking into account the effect of protein structure flexibility,³⁰ we could not consider such a dynamic property in the virtual screening of DUSP16 inhibitors because the receptor had to be treated as a rigid body in the current version of the AutoDock program. Of the 260,000 compounds screened with docking simulations in the active site of the homology-modeled DUSP16 structure, 150 top-scored compounds were selected as the virtual hits. A total of 148 of them were commercially available from the compound supplier and tested for inhibitory activity against DUSP16 by in vitro enzyme assay. As a result, we found seven compounds that inhibited the catalytic activity of DUSP16 by more than 50% at the concentration of 20 µM, which were selected to measure the IC₅₀ values. The molecular structures and the inhibitory activities of the newly found inhibitors are shown in Figure 2 and Table 1 , respectively. The dose-response behaviors of the inhibitors measured to obtain their IC₅₀ values are presented in Supplementary Figure S2. We note that compounds 1–7 exhibit good inhibitory activities, with the associated IC₅₀ values ranging from 1.3 to 21.3 µM. Benzoate and nitrobenzene moieties included in the molecular structures of 1–7 are expected to serve as a surrogate for the phosphotyrosine group of the protein substrates. It is also a common structural feature for 1–7 to include multiple nonpolar aromatic groups. Therefore, the newly identified inhibitors seem to be capable of establishing strong hydrogen bonds and van der Waals contact as well in the active site of DUSP16. Because 1–7 were also computationally screened for having desirable physicochemical properties as a drug candidate, they deserve consideration for further development by structure-activity relationship (SAR) methods to optimize the inhibitory activity against DUSP16.

Figure 2.

Chemical structures of seven DUSP16 inhibitors identified with virtual screening.

Table 1.

IC₅₀ values (in µM) of 1–7 with respect to DUSP16, DUSP10, and DUSP25.^a

Compound	DUSP16	DUSP10	DUSP25
1	1.3 (±0.2)	4.2 (±0.6)	4.3 (±1.1)
2	4.4 (±0.7)	8.4 (±2.0)	7.4 (±0.8)
3	4.8 (±1.2)	1.7 (±0.5)	0.2 (±0.1)
4	7.2 (±2.9)	13.7 (±0.8)	16.0 (±5.9)
5	10.9 (±4.7)	36.4 (±7.1)	34.2 (±6.6)
6	14.8 (±3.8)	27.2 (±6.3)	19.2 (±4.9)
7	21.3 (±5.8)	28.0 (±1.4)	8.5 (±1.7)

Each IC₅₀ value was obtained by averaging the two data points determined from the dose-response behaviors of the inhibitor measured in duplicate. Standard deviations are presented in parentheses.

Because 1–7 included the delocalized electronic systems with multiple conjugated aromatic moieties, we addressed the possibility that the inhibitors would fluorescently interfere with the assay signal. This validation study started with the incubation of the enzymatic reaction mixture without the treatment of inhibitor to obtain the maximum fluorescence signal. At this point, the enzymatic reaction was terminated by the treatment of vanadate, which was followed by the addition of each inhibitor in the reaction mixture to determine its effect of fluorescence interference. The extent of decrease in fluorescence signal was then measured at the inhibitor concentration of 5 µM. The fluorescence signals appeared to be reduced by only 0% to 3% because of the interference effects from the inhibitors. These small interference effects indicate that the major contribution to the decrease in fluorescence signal in enzyme assays should come from the inhibition of enzymatic activity rather than from the role of the fluorescence quencher played by the inhibitors.

Because the selectivity has been one of the most important issues in the development of phosphatase inhibitors, we determined the inhibitory activities of 1–7 with respect to the two other DUSPs, DUSP10 and DUSP25, because they could be purified in the active form to the amount required in the enzyme inhibition assay. As shown in Table 1 and in Supplementary Figure S3, most newly identified inhibitors reveal lower potencies for DUSP10 and DUSP25 than for DUSP16, although 3 appears to be also a potent inhibitor of the former. Thus, the inhibitors 1, 2, and 4–7 seem to bind in the active site of DUSP16 in a stronger fashion than in those of the other structurally similar DUSPs. The significant differences in the inhibitory activities for the three DUSPs also indicate that 1–7 would impair their phosphatase activity through the nonbond interactions in the active site instead of making a covalent bond with the side-chain thiolate ion of the catalytic cysteine residue contained in common.

Related to the imperfections in the simple desolvation model implemented in the original AutoDock program, we found that the rankings of 2, 4, and 5 went beyond 150 in the virtual screening with the previous scoring function. Therefore, only four inhibitors would be found in virtual screening without the modification of the desolvation term, whereas all seven inhibitors 1–7 could be identified as a DUSP16 inhibitor with the improved scoring function. The increase in the number of the actual inhibitors from four to seven can thus exemplify the superiority of the modified scoring function over the previous one.

To gain structural insight into the inhibition mechanisms for the identified DUSP16 inhibitors, we examined their binding modes in the active site in a comparative fashion. Figure 3 displays the lowest-energy conformations of 1–7 in the active site of DUSP16 calculated with the modified AutoDock program. The results of these docking simulations are self-consistent because the functional groups of similar chemical character are placed in similar ways with comparable interactions with the protein groups. As revealed by the superposition of the docked structures, for example, the polar moieties of 1–7 that serve as a surrogate for the substrate phosphotyrosine group point toward the PTP loop including the catalytic cysteine residue (Cys244), whereas their hydrophobic groups are directed to the Q loop connecting the two alpha helices near the active site. These common features in the calculated binding modes indicate that a potent DUSP16 inhibitor should include an effective surrogate for the substrate phosphortyrosine group and simultaneously the hydrophobic moieties for binding to the Q loop. To address the possibility of the allosteric inhibition of DUSP16 by the identified inhibitors, docking simulations were carried out with the grid maps for the receptor model so as to include the entire phosphatase domain. However, the binding configuration in which an inhibitor resides outside the active site was not observed for any of the inhibitors found in this study. These results support the possibility that the inhibitors would impair the catalytic activity of DUSP16 through binding in the active site rather than in the peripheral binding sites.

Figure 3.

Comparative view of the binding modes of 1–7 in the active site of DUSP16. Carbon atoms of 1–7 are indicated in green, cyan, black, gray, pink, orange, and violet, respectively. The positions of PTP loop, WPD and Q loops, and Cys244 are also indicated.

The calculated binding mode of 1 in the active site of DUSP16 is illustrated in Figure 4 . The inhibitor is found to be in close contact with Cys244–Arg250 and Ile281–Phe287, which belong to the PTP and Q loops, respectively. We note that one of the terminal carboxylate oxygens of 1 plays the role of hydrogen bond receptor with respect to the side-chain guanidinium ion of Arg250 and the backbone amidic groups of Cys244, Ala246, and Gly247 at the bottom of the active site. The other carboxylate oxygen also establishes a hydrogen bond with backbone amidic nitrogen of Ser249. It is apparent that these five hydrogen bonds involving the carboxylate ion can serve as an anchor for positioning the inhibitor in the active site. It is also noteworthy that the carboxylate carbon of 1 resides in proximity to the side-chain thiolate ion of Cys244 with the associated interatomic distance of 4.06 Å. Judging from the proximity to Cys244 and the establishment of multiple hydrogen bonds in the active site, the benzoate moiety of 1 seems to play the role of an effective surrogate for the substrate phosphotyrosine group. An additional hydrogen bond is found between the aminocarbonyl oxygen on the pyrrol-2-one ring of 1 and the side-chain hydroxyl group of Ser282, which should also be a significant binding force in the DUSP16-1 complex. The inhibitor 1 appears to be further stabilized in the active site of DUSP16 by the hydrophobic interactions of its nonpolor groups with the side chains of Ala246, Ile248, Pro283, and Phe287. Thus, the overall structural features of the DUSP16-1 complex obtained from docking simulations indicate that the micromolar inhibitory activity of 1 would stem from the multiple hydrogen bonds and hydrophobic interactions established simultaneously in the active site.

Figure 4.

Calculated binding mode of 1 in the active site of DUSP16. Carbon atoms of the protein and the ligand are indicated in cyan and green, respectively. Each dotted line indicates a hydrogen bond.

Figure 5 shows the lowest-energy binding mode of 2 in the active site of DUSP16. The binding mode of 2 is similar to that of 1 in that the roles of hydrogen bond donor with respect to the terminal 2-nitrophenol moiety of 2 are played by the side-chain guanidinium ion of Arg250 and the backbone amidic groups located at the bottom of the active site. Like the carboxylate group of 1, the nitrogen atom of the terminal nitro group of 2 resides in the vicinity of the side-chain thiolate ion of Cys244 at the distance of 4.02 Å. Therefore, the 2-nitrophenol group of 2 is expected to serve as an effective surrogate for the substrate phosphotyrosine group. A stable hydrogen bond is also established between the nitrile moiety of 2 and side-chain hydroxyl group of Ser249, which may also contribute to the stabilization of the inhibitor in the active site. Hydrophobic interactions in the DUSP16-2 complex appear to be established in a similar fashion to those in the DUSP16-1 complex: the nonpolar aromatic groups of 2 form the van der Waals contacts with the side chains of Ala246, Ile248, Pro283, and Phe287. Although the inhibitory activity of 2 is a little lower than that of 1 ( Table 1 ), 2 is expected to serve as a good inhibitor scaffold from which much more potent inhibitors can be derivatized because of its relatively low molecular weight of ~393.

Figure 5.

Calculated binding mode of 2 in the active site of DUSP16. Carbon atoms of the protein and the ligand are indicated in cyan and green, respectively. Each dotted line indicates a hydrogen bond.

As can be inferred from the overlaid docking poses in Figure 3 , 3–7 exhibit a binding mode similar to those of 1 and 2 in that the terminal carboxylate/nitro moieties and the aromatic rings establish the hydrogen bonds with the active-site residues and the van der Waals contacts with nonpolar residues in PTP, Q, and WPD loops, respectively. Thus, it seems to be necessary for the potent DUSP16 inhibitors to have a proper polar group at the end of the chemical structure to form the multiple hydrogen bonds in the active site and simultaneously the nonpolar groups to be stabilized through the hydrophobic interactions with nonpolar resides in the loop structures.

In conclusion, we have identified seven novel inhibitors of DUSP16 by means of a computer-aided drug design protocol involving the homology modeling and the structure-based virtual screening. To increase the accuracy in the calculated binding affinity between DUSP16 and the putative inhibitors, the scoring function was improved in such a way as to reflect the effect of ligand solvation on protein-ligand docking. The newly found inhibitors revealed a good potency, with IC₅₀ values ranging from 1.3 to 21.3 µM. Because these inhibitors were also found to have significant inhibitory activity for other DUSPs, future SAR studies to develop the lead compounds for a new anticancer medicine should be carried out in such a way as to improve the selectivity in the inhibition of DUSP16 as well as to optimize the inhibitor potency. Detailed binding mode analyses with docking simulations indicated that the inhibitors could be stabilized through the multiple hydrogen bonds and the van der Waals contacts established simultaneously in the active site of DUSP16.

Footnotes

Acknowledgements

The authors thank the Olson group at Scripps Research Institute for kindly providing the AutoDock program.

Supplementary material for this article is available on the Journal of Biomolecular Screening Web site at .

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2011-0022858), and a Hanyang University internal grant.

References

Chang

Karin

Mammalian MAP Kinase Signalling Cascades. Nature 2001, 410, 37–40.

Patterson

K. I.

Brummer

O’Brien

P. M.

. Dual-Specificity Phosphatases: Critical Regulators with Diverse Cellular Targets. Biochem. J. 2009, 418, 475–489.

Musikacharoen

Bandow

Kakimoto

. Functional Involvement of Dual Specificity Phosphatase 16 (DUSP16), a c-Jun N-Terminal Kinase-Specific Phosphatase, in the Regulation of T Helper Cell Differentiation. J. Biol. Chem. 2011, 286, 24896–24905.

Lee

Syed

Taylor

. DUSP16 Is an Epigenetically Regulated Determinant of JNK Signalling in Burkitt’s Lymphoma. Br. J. Cancer 2010, 103, 265–274.

Zhu

Shen

Gao

. PI3K/Akt-Independent Negative Regulation of JNK Signaling by MKP-7 after Cerebral Ischemia in Rat Hippocampus. BMC Neurosci. 2013, 14, 1.

Nizamutdinova

I. T.

Kim

Y. M.

Lee

J. H.

. MKP-7, a Negative Regulator of JNK, Regulates VCAM-1 Expression through IRF-1. Cell. Signal. 2012, 24, 866–872.

Cybulsky

M. I.

Iiyama

. A Major Role for VCAM-1, but Not ICAM-1, in Early Atherosclerosis. J. Clin. Invest. 2001, 107, 1255–1262.

Almo

S. C.

Bonanno

J. B.

Sauder

J. M.

. Structural Genomics of Protein Phosphatases. J. Struct. Funct. Genom. 2007, 8, 121–140.

Lavecchia

Di Giovanni

Virtual Screening Strategies in Drug Discovery: A Critical Review. Curr. Med. Chem. 2013, 20, 2839–2860.

10.

Lill

Virtual Screening in Drug Design. Methods Mol. Biol. 2013, 993, 1–12.

11.

Huang

S. Y.

Grinter

S. Z.

Zou

Scoring Functions and Their Evaluation Methods for Protein-Ligand Docking: Recent Advances and Future Directions. Phys. Chem. Chem. Phys. 2010, 12, 12899–12908.

12.

Shoichet

B. K.

Leach

A. R.

Kuntz

I. D.

Ligand Solvation in Molecular Docking. Proteins 1999, 34, 4–16.

13.

Thompson

J. D.

Higgins

D. G.

Gibson

T. J.

CLUSTAL W: Improving the Sensitivity of Progressive Multiple Sequence Alignment through Sequence Weighting, Position-Specific Gap Penalties and Weight Matrix Choice. Nucleic Acids Res. 1994, 22, 4673–4680.

14.

Sali

Blundell

T. L.

Comparative Protein Modelling by Satisfaction of Spatial Restraints. J. Mol. Biol. 1993, 234, 779–815.

15.

Fiser

R. K.

Sali

Modeling of Loops in Protein Structures. Protein Sci. 2000, 9, 1753–1773.

16.

Sippl

M. J.

Recognition of Errors in Three-Dimensional Structures of Proteins. Proteins 1993, 17, 355–362.

17.

Jeffrey

G. A.

An Introduction to Hydrogen Bonding; Oxford University Press: Oxford, 1997.

18.

Lipinski

C. A.

Lombardo

Dominy

B. W.

. Experimental and Computational Approaches to Estimate Solubility and Permeability in Drug Discovery and Development Settings. Adv. Drug Deliv. Rev. 1997, 23, 3–20.

19.

Gasteiger

Marsili

Iterative Partial Equalization of Orbital Electronegativity: A Rapid Access to Atomic Charges. Tetrahedron 1980, 36, 3219–3228.

20.

Morris

G. M.

Goodsell

D. S.

Halliday

R. S.

. Automated Docking Using a Lamarckian Genetic Algorithm and an Empirical Binding Free Energy Function. J. Comput. Chem. 1998, 19, 1639–1662.

21.

Mehler

E. L.

Solmajer

Electrostatic Effects in Proteins: Comparison of Dielectric and Charge Models. Protein Eng. 1991, 4, 903–910.

22.

Stouten

P. F. W.

Frömmel

Nakamura

Sander

An Effective Solvation Term Based on Atomic Occupancies for Use in Protein Simulations. Mol. Simul. 1993, 10, 97–120.

23.

Choi

Kang

Park

New Solvation Free Energy Function Comprising Intermolecular Solvation and Intramolecular Self-Solvation Terms. J. Cheminform. 2013, 5, 8.

24.

Park

Hong

Kim

. Discovery of Picomolar ABL Kinase Inhibitors Equipotent for Wild Type and T315I Mutant via Structure-Based De Novo Design. J. Am. Chem. Soc. 2013, 135, 8227–8237.

25.

Park

K. R.

. Identification of Novel PTPRQ Phosphatase Inhibitors Based on the Virtual Screening with Docking Simulations. Theor. Biol. Med. Model. 2013, 10, 49.

26.

Welte

Baringhaus

K.-H.

Schmider

. 6,8-Difluoro-4-Methylumbiliferyl Phosphate: A Fluorogenic Substrate for Protein Tyrosine Phosphatases. Anal. Biochem. 2005, 338, 32–38.

27.

Park

Kyung

Lee

H.-J.

. Virtual Screening and Biochemical Evaluation of the Inhibitors of Dual-Specificity Phosphatase 26. Med. Chem. Res. 2013, 22, 3905–3910.

28.

Park

Bahn

Y. J.

Jung

S.-K.

. Discovery of Novel Cdc25 Phosphatase Inhibitors with Micromolar Activity Based on the Structure-Based Virtual Screening. J. Med. Chem. 2008, 51, 5533–5541.

29.

Baker

Sali

Protein Structure Prediction and Structural Genomics. Science 2001, 294, 93–96.

30.

Nichols

S. E.

Baron

Ivetac

. Predictive Power of Molecular Dynamics Receptor Structures in Virtual Screening. J. Chem. Inf. Model. 2011, 51, 1439–1446.