Discovery of Novel DUSP4 Inhibitors through the Virtual Screening with Docking Simulations

Abstract

Dual specificity protein phosphatase 4 (DUSP4) has been considered a promising target for the development of therapeutics for various human cancers. Here, we report the first example for a successful application of the structure-based virtual screening to identify the novel small-molecule DUSP4 inhibitors. As a consequence of the virtual screening with the modified scoring function to include an effective molecular solvation free energy term, five micromolar DUSP4 inhibitors are found with the associated $IC_{50}$ values ranging from 3.5 to $10.8{\mu}M$. Because these newly identified inhibitors were also screened for having desirable physicochemical properties as a drug candidate, they may serve as a starting point of the structure-activity relationship study to optimize the medical efficacy. Structural features relevant to the stabilization of the new inhibitors in the active site of DUSP4 are discussed in detail.

Introduction

Dual-specificity protein phosphatases (DUSPs) hydrolyze the phosphorylated serine/threonine and tyrosine residues, and thereby play vital roles in immune activation, brain function, and cell-growth signalling. The structural superposition of all human DUSPs exhibits a good alignment in the core region with some deviations in the loops and secondary structures outside the core. Among a variety of human DUSPs, DUSP4 (also known as MKP2) dephosphorylates various MAP kinases including ERK1/2, p38, and JNK. DUSP4 expression was found to be up-regulated in human melanoma cell lines1 and ovarian cancer.2 This is consistent with the earlier finding that deletion of DUSP4 by RNAi caused the proliferative inhibition of a mammary carcinoma cell line.3 Recently, Gröschl et al. reported that the overexpression of DUSP4 was also associated with microsatellite instability in colorectal cancer and observed in most colorectal cancer patients,4 which indicated its significant role in the promotion of cancer. The carcinoma developed in the RET-MEN2A transgenic mice that carried the RET kinase gene with multiple endocrine neoplasia 2A mutation also showed the increased expression of DUSP4.5 The involvement of DUSP4 in caner was further evidenced by the experimental finding that the knockout of DUSP4 gene in mouse led to the significant retardation of cell proliferation and to the increase in apoptosis.6 Thus, the inhibition of DUSP4 activity can be a promising strategy for the development of anti-cancer medicines.

Three dimensional (3D) structure of the catalytic domain of DUSP4 was reported in complex with a substrate analog (sulfate ion).7 The overall folding of DUSP4 is similar to those of the typical dual-specificity phosphatases with a central twisted five stranded β-sheets surrounded by six α-helices. The residues in the PTP loop (C280-Q281-A282-G283-I284-S285) and the neighboring general acid loop constitute the active site pocket in a flat and shallow form, which is consistent with the broad substrate specificity of DUSP4. Structural information about the nature of the active site and the interactions with a small-molecule ligand can be very helpful in designing the potent inhibitors that may develop into a drug candidate. Indeed, the usefulness of such structural information has been well appreciated in designing the potent and selective small-molecule inhibitors of various phosphatases.8 Nonetheless, the discovery of DUSP4 inhibitors has lagged behind the biological and structural studies. To the best of our knowledge, no small-molecule inhibitor has been in reported so far at least in the literature.

In the present study, we aim to identify the potent DUSP4 inhibitors by means of a structure-based drug design protocol involving the virtual screening with docking simulations and in vitro enzyme assay. The characteristic feature that discriminates our virtual screening approach from the others lies in the implementation of a molecular solvation free energy term in the scoring function to calculate the binding free energy between DUSP4 and the putative ligands, which would have an effect of increasing the probability of finding the actual inhibitors in enzyme assays.9 We find in this study that the docking simulations with the improved binding free energy function can be a useful computational tool for enriching the chemical library with molecules that are likely to have inhibitory activity against DUSP4 as well as for elucidating the potencies of the identified inhibitors.

 

Computational and Experimental Methods

The X-ray crystal structure of human DUSP4 in complex with the sulfate ion used as a substrate analogue (PDB entry: 3EZZ)7 was selected as the receptor model in the virtual screening with docking simulations. After removing the solvent molecules, hydrogen atoms were added to each protein atom. A special attention was paid to assign the protonation states of the ionizable Asp, Glu, His, and Lys residues in the original X-ray crystal structure of DUSP4. The side chains of Asp and Glu residues were assumed to be neutral if one of their carboxylate oxygens pointed toward a hydrogen-bond accepting group such as the backbone aminocarbonyl oxygen at a distance within 3.5 Å, a generally accepted distance limit for a hydrogen bond of moderate strength.10 Similarly, the lysine side chains were assumed to be protonated unless the NZ atom was in proximity to a hydrogen-bond donating group. The same procedure was also applied to determine the protonation states of ND and NE atoms in His residues.

The docking library for DUSP4 comprising about 260,000 compounds was constructed from the latest version of the chemical database distributed by Interbioscreen (http://www.ibscreen.com) containing approximately 480,000 synthetic and natural compounds. Prior to the virtual screening with docking simulations, they were filtrated on the basis of Lipinski’s “Rule of Five” to adopt only the compounds with the physicochemical properties of potential drug candidates11 and without reactive functional group(s). All of the compounds included in the docking library were then processed with the CORINA program to generate their 3D atomic coordinates, followed by the assignment of Gasteiger-Marsilli atomic charges.12 We used the AutoDock program13 in the virtual screening because the outperformance of its scoring function over those of the others had been shown for various target proteins.14-17 AMBER force field parameters were assigned to compute the van der Waals interactions and the internal energy of a ligand. Docking simulations with AutoDock were then carried out in the active site of DUSP4 to score and rank the compounds in the docking library according to their calculated binding affinities.

In the actual docking simulations between DUSP4 and the putative ligands, we used the empirical AutoDock scoring function improved by the implementation of a proper molecular solvation free energy function. The modified scoring function can be written as follows.

Here WvdW, Whbond, Welec, Wtor, and Wsol are the weighting factors of van der Waals, hydrogen bond, electrostatic interactions, torsional term, and desolvation energy of inhibitors, respectively. rij represents the interatomic distance, and Aij, Bij, Cij, and Dij are related to the depths of the potential energy well and the equilibrium separations between the protein and ligand atoms. The hydrogen bond term has an additional weighting factor, E(t), representing the angledependent directionality. Cubic equation approach was applied to obtain the dielectric constant required in computing the interatomic electrostatic interactions between DUSP4 and a ligand molecule.18 In the entropic term, Ntor is the number of rotatable bonds in the ligand. In the desolvation term, Si and Vi are the solvation parameter and the fragmental volume of atom i,19 respectively, while Occimax stands for the maximum atomic occupancy. In the calculation of the 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.20 The accuracy of virtual screening seems to be enhanced due to the inclusion of this solvation free energy term in the scoring function because the underestimation of ligand solvation effects often leads to the overestimation of the binding affinity of a ligand with many polar atoms.9 Indeed, the superiority of this modified scoring function to the previous one was demonstrated in recent studies for virtual screening of kinase and phosphatase inhibitors.21-24

The catalytic domain of DUSP4 was produced by using an E. coli over-expression system. Briefly, the pET28a vector containing the DUSP4 catalytic domain was transformed to BL21 competent cells. Then the cells were induced with 0.1 mM IPTG at 18 ℃ overnight. Cell pellets were harvested and lysed with a buffer containing 50 mM Tris-Cl (pH 7.5), 500 mM NaCl, 1 mM phenylmethylsulfonyl fluoride (PMSF), 0.05％ 2-mercaptoethanol, and 5％ glycerol. The DUSP4 catalytic domain was purified by Ni-NTA affinity and Qsepharose chromatographies. The purified protein was dialyzed against 20 mM Tris-Cl (pH 8.0), 50 mM NaCl, 10 mM 2-mercaptoethanol, 0.5 mM EDTA, and 5％ glycerol.

A total of 143 compounds were tested for inhibitors of DUSP4. The phosphatase assays were performed by using 6,8-difluoro-4-methyl-umbelliferyl phosphate (DiFMUP) as a fluorogenic substrate. The concentrations of DUSP4 (500 nM) and DiFMUP (10 µM) was determined by selecting a concentration showing the typical Michaelis-Menten kinetics. The purified DUSP4 (500 nM) and a candidate inhibitor (20 µM) were mixed in a reaction buffer containing 20 mM Tris-Cl (pH 8.0), 0.01％ Triton X-100, and 5 mM DTT. The reaction was continued for 20 min at room temperature. The enzymatic activity was measured using the Perkin Elmer 2030 instrument with a wavelength of 355 and 460 nm as excitation and emission wavelengths, respectively. Initial velocities of reactions were estimated for various concentrations (0-50 µM) of inhibitors. The inhibitory activities were measured in duplicate at the concentrations of 0, 0.625, 1.25, 2.5, 5, 10, 20, and 50 µM. The IC50 value of each inhibitor was measured from the direct regression analysis by using a four parameter sigmoidal curve in the OriginPro program.

 

Results and Discussion

Of the 260,000 compounds screened with docking simulations, 150 top-scored compounds were selected as virtual hits. 143 of them were available from the compound supplier and were tested for inhibitory activity against DUSP4 by in vitro enzyme assays. As a result, we identified the five compounds that inhibited the catalytic activity of DUSP4 by more than 50％ at the concentration of 10 µM, which were selected to determine the IC50 values. The chemical structures and IC50 values of the newly identified inhibitors are shown in Figure 1. We note that compounds 1-5 have good potency with the associated IC50 values ranging from 3.5 to 10.8 µM. Benzoate moiety in 1 and 5, nitrobenzene in 2 and 3, and sulfonylurea group in 4 seem to serve as an effective surrogate for the substrate phosphotyrosine group. These polar groups are therefore expected to be capable of establishing the strong multiple hydrogen bonds with the amino acid residues in the active site. Besides the good inhibitory activities, 1-5 were also screened for having desirable physicochemical properties as a drug candidate. Therefore, they deserve consideration for further development by the structure-activity relationship (SAR) studies to optimize the anticancer activity.

Figure 1.Chemical structures and IC50 values of the newly identified DUSP4 inhibitors.

To estimate the accuracy of the scoring function used in virtual screening, we calculated the Pearson coefficient (R) between the calculated binding free energies of 1-5 and their corresponding experimental IC50 values. The calculated R value amounted to 0.84, which indicated that the experimental and computational results would be closely related.

To further assess the merit of the modified scoring function implementing the molecular solvation term, we compared the binding free energies of 260,000 compounds calculated with the modified scoring function with those obtained with the previous one that lacks the desolvation term. The results show that the rankings of 2 and 5 go beyond 150 when the old previous scoring function was employed to estimate the binding affinities. This indicates that only 1, 3, and 4 could be found in virtual screening with the previous scoring function, and that all of the five inhibitors could be identified as DUSP4 inhibitor only when the solvation term was included in the scoring function. This result exemplifies the superiority of the scoring function with the desolvation term to the previous one.

It is worth noting that the selectivity is one of the most important properties DUSP inhibitors should have as a drug candidate because of the high structural similarities of the active sites among various DUSPs. Therefore, we also measured the inhibitory activities of 1-5 for DUSP16, which is one of the highly homologous proteins with respect to DUSP4. These inhibition assays for selectivity were done in duplicates at the concentration of 10 µM. The results show that 1, 2, and 3 exhibit also a significant inhibitory activity against DUSP16, which exemplifies the difficulty in the discovery of specific DUSP inhibitors. The simultaneous inhibitions of DUSP4 and DUSP16 by 1-3 are actually not surprising because they share a highly conserved catalytic module. To obtain the specific inhibitors for DUSP4, therefore, it seems that some chemical groups should be added to 1-3 in such a way that the resulting derivatives can be stabilized not only in the active site but also in other peripheral binding pockets.

To obtain structural insight into the inhibitory mechanisms of the identified DUSP4 inhibitors, their binding modes in the active site were investigated in a comparative fashion. Figure 2 shows the lowest-energy conformations of 1-5 in the active site gorge of DUSP4 calculated with the modified scoring function. The results of docking simulations are selfconsistent 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 hydrophilic moieties that serve as a surrogate for the substrate phosphotyrosine group are bound in the PTP loop including the catalytic cysteine residue (Cys280) while the hydrophobic groups point toward the loops located above the active site. This common feature in the binding modes indicates the necessity for an effective surrogate of the substrate phosphotyrosine group and simultaneously for the hydrophobic groups for binding to the loop structures as the key structural components of potent DUSP4 inhibitors. In order to examine the possibility of the allosteric inhibition of DUSP4 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 of DUSP4. However, the binding configuration in which an inhibitor resides outside the active site was not observed for any of the new inhibitors. These results support the possibility that the inhibitors would impair the catalytic activity of DUSP4 through the specific binding in the active site.

Figure 2.Comparative view of the binding modes of 1-5 in the active site of DUSP4. Carbon atoms of 1-5 are indicated in green, cyan, black, pink, and gray, respectively. The positions of PTP loop, WPD loop, and Cys280 are also indicated.

The calculated binding mode of 1 in the active site of DUSP4 is shown in Figure 3. We note that two carboxylate oxygens of 1 play the role of hydrogen bond receptor with respect to the backbone amidic nitrogens of Ala282, Ile284, Ser285, and Arg286 at the bottom of active site. Apparently, these four hydrogen bonds seem to play a critical role in anchoring the inhibitor in the active site. It is also noted that the carboxylate carbon of 1 resides in the vicinity of the sidechain thiolate ion of Cys280 with the associated interatomic distance of 3.8 Å. Judging from the proximity to Cys280 and the formation of multiple hydrogen bonds at the active site, the benzoate moiety of 1 should be an effective surrogate for the phosphotyrosine group in the substrates of DUSP4. The carbonyl oxygen on the central lactone ring of 1 appears to form the additional hydrogen bonds with the side-chain guanidinium ion of Arg286 and the backbone amide moiety of Ser230 in the bifurcated form, which should also be a significant binding force in the DUSP4-1 complex. The inhibitor 1 can be further stabilized in the active site of DUSP4 by establishing the hydrophobic interactions between its nonpolor groups and the side chains of Ala208, Cys232, Pro233, Ala282, and Ile284. Thus, the overall structural features derived from docking simulations indicate that the micromolar inhibitory activity of 1 should stem from the combined effect of the multiple hydrogen bonds and hydrophobic interactions established simultaneously in the active site of DUSP4.

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

Figure 4 shows the lowest-energy binding mode of 2 in the active site of DUSP4. In this calculated structure of DUSP4-2 complex, the roles of hydrogen bond donor with respect to the terminal nitrobenzene group of 2 are played by the backbone amidic groups of Ile284 and Ser285. Similar to the benzoate moiety in the DUSP4-1 complex, the terminal nitro group of 2 stays close to the side-chain thiolate ion of Cys280 at the interatomic distances within 4.1 Å. This indicates that it can also serve as an effective surrogate for the substrate phosphotyrosine group that is necessary for the effective inhibition of DUSP4. Two additional hydrogen bonds are formed in the DUSP4-2 complex: one between the central aminocarbonyl oxygen of 2 and the side-chain guanidinium ion of Arg286, and the other between the imine group of 2 and the side-chain carboxylate ion of Asp231. Because these two hydrogen bonds involve a charged group, they seem to play more significant role in stabilizing the inhibitor in the active site than those established between neutral groups. Hydrophobic interactions in DUSP4-2 complex appear to be established in the similar form to those in DUSP-1 complex: its phenyl, furan, and pyridine rings form a van der Waals contact with the side chains of Ala208, Cys232, Pro233, Ala282, and Ile284. However, the hydrogenbond interactions in DUSP4-2 complex seem to be weaker than those in DUSP4-1 one because the number of hydrogen bonds decreases from six in the former to four in the latter. This can be invoked to explain a little lower inhibitory activity of 2 than 1.

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

Binding modes of 3-5 are similar to that of 1 and 2 in that the terminal polar groups form the multiple hydrogen bonds with the backbone amidic groups on the PTP loop in the vicinity of Cys280 at the bottom of active site. This confirms the usefulness of the benzoate, nitrobenzene, and sulfonylurea groups as a surrogate for the substrate phosphotyrosine group. As can be inferred from the overlaid docked structures of the inhibitors in Figure 2, however, the hydrophobic interactions in DUSP4-3 and DUSP4-4 complexes are established in a different form from those in DUSP4-1, DUSP4-2, and DUSP4-5 complexes. This indicates the presence of two hydrophobic binding pockets in which nonpolar groups of the inhibitors can be stabilized. Therefore, the inhibitory activities of 1-5 seem to be enhanced by the addition of chemical groups in such a way that the resulting derivatives can bind in both hydrophobic binding pockets.

 

Conclusions

In summary, we have identified five novel inhibitors of DUSP4 by applying a computer-aided drug design protocol involving the structure-based virtual screening with docking simulations under consideration of the effects of ligand solvation on the protein-ligand binding affinity. These inhibitors revealed a high potency with the IC50 values ranging from 3.5 to 10.8 µM and were also screened for having desirable physicochemical properties as a drug candidate. Therefore, each of the newly identified inhibitors deserves consideration for further development by SAR studies to optimize the anticancer activity. Detailed binding mode analyses with docking simulations showed that the inhibitors could be stabilized in active site by the simultaneous establishment of multiple hydrogen bonds and van der Waals contacts.