Transcription factor STAT3 is an important target protein that is involved in a multitude of human cancers. In this work, we focused on a specific set of 12 peptidomimetic compounds that mimic the pTyr-Xaa-Yaa-Gln recognition motif and were designed to bind with the SH2 domain of STAT3 and prevent its dimerization which is a critical event leading up to the transcription of anti-apoptotic genes. Experimental binding affinities of the peptidomimetics were measured using fluorescence polarization and a range of affinity values were observed for the 12 peptidomimetics. Binding affinities for the peptidomimetics, expressed as IC50 values, range from 39 nM for a strong binder to over 100,000 nM for a weak binder. Since experimental structures of the complexes formed between the peptidomimetics and the SH2 domain are unavailable, we used a computational strategy to model the complexes.
Our modeling strategy proceeded in two steps. In the first, we generated docked conformations of the peptidomimetics using a computational AutoDock-based incremental docking protocol that was developed by us for docking large compounds in a fast and accurate manner 
. The peptidomimetics in our dataset are all large compounds with the number of rotatable bonds ranging from 9 to 22. In the second step of our modeling strategy, we selected the best docked conformation and then ran molecular dynamics simulations of the complex in a solvated box. Molecular dynamics simulations served multiple purposes. The flexibility of the SH2 domain was taken into account, fluctuations of the bound conformations over the length of molecular dynamics simulation were computed, and finally, rigorous estimates of binding affinities, as a sum of normal-mode analysis based entropic component and MMPB/GBSA based non-entropic component, were computed. Accurate estimates of binding affinities are very important for differentiating strong binders from weak binders, and therefore, a positive correlation between the experimental binding affinities and estimated binding affinities is desired. Our two-step modeling strategy resulted in a high positive correlation (R
0.63) between the experimental and estimated affinities.
For each of the 12 peptidomimetics, we performed molecular dynamics simulations for a production length of 10 ns. The trajectory data for each simulation was output at 10 ps. Thus, we obtained 1000 conformations for each peptidomimetic in complex with the SH2 domain. The average fluctuation of the conformations of each peptidomimetic was measured as RMSF (root mean square fluctuation) values. The weak binders displayed larger fluctuation as compared to the strong binders. A clustering of the conformations showed the preferred binding modes of the peptidomimetics. Three strong binders, with IC50 values equal to 190 nM (comp70), 83 nM (comp134), and 68 nM (comp121), displayed three different but stable binding modes: the bent mode, the extended mode, and the wedged mode respectively. The peptidomimetics in these three binding modes showed very small (<1.0 Å) conformational fluctuations in the molecular dynamics simulations, a large number of stable hydrogen bond interactions with the residues of the SH2 domain, and the estimated binding affinities value were low in accordance with the experimental binding affinities.
Previous modeling studies related to SH2 domain binding have proposed the bent and the extended binding modes 
. In this paper, we propose a new binding mode which we term the wedged mode. In the wedged mode, the peptidomimetic (comp121) binds to the SH2 domain such that the negatively charged phosphate group of the pTyr residue sits inside a pocket which has a positive electrostatic potential and the C-terminal benzyl group gets wedged in a cavity formed by two loops of the SH2 domain described by the residues 623–629 and 656–668 respectively. Apart from the stable hydrogen bond interactions with the residues in the phosphate-binding pocket, hydrogen bonds also exist between the peptidomimetic and residues on the two loops. The RMSF value for the 1000 conformations of the comp121 is 0.91 Å and is the lowest among the RMSF values for the 12 peptidomimetics.
Despite the overall success of modeling strategy as described in this paper, there were exceptions to the observed trends. For example, in the case of comp140 which is a relatively strong binder (IC50
105 nM), we obtained a large RMSF value and estimated binding affinities that are comparable to those of weak binders. This anomaly could be attributed to an inaccurate starting docked conformation of the peptidomimetic. In the molecular dynamics simulation, an inaccurate starting docked conformation would result in trajectory that leads to inaccurate estimation of binding affinity. It should be noted that computational docking of large ligands such as peptidomimetics in our dataset is extremely challenging. Although our incremental docking protocol improves docking of large ligands, more work needs to be done in this area.
The computational modeling strategy described in this paper and the subsequent data analysis, nonetheless, reveals important aspects of the peptidomimetic binding to the SH2 domain of STAT3. Not only were we able to estimate binding affinities that were well correlated with experimental binding affinities, we were also able to identify binding modes, including a novel wedge mode, that result in stable binding interactions. A typical peptidomimetic drug design process that is based on a specific motif involves designing peptidomimetics with diverse chemical modifications. Accurate estimation of binding affinities using our method could help in predicting which modifications could lead to strong binding. The knowledge gained by this study could also be used to improve the design of the peptidomimetics by better targeting the sub-binding-pockets identified in this paper with structural modifications or conformational restraints. The proposed novel wedge binding mode could prove very useful in this regard.