If You Cannot Win Them, Join Them: Understanding New Ways to Target STAT3 by Small Molecules.

Signal transducer activator of transcription 3 (STAT3) is among the most investigated oncogenic transcription factors, as it is highly associated with cancer initiation, progression, metastasis, chemoresistance, and immune evasion. Evidences from both preclinical and clinical studies have demonstrated that STAT3 plays a critical role in several malignancies associated with poor prognosis such as glioblastoma and triple-negative breast cancer, and STAT3 inhibitors have shown efficacy in inhibiting cancer growth and metastasis. Constitutive activation of STAT3 by mutations occurs frequently in tumor cells and directly contributes to many malignant phenotypes. Unfortunately, detailed structural biology studies on STAT3 as well as target-based drug discovery efforts have been hampered by difficulties in the expression and purification of the full-length STAT3 and a lack of ligand-bound crystal structures. Considering these, molecular modeling and simulations offer an attractive strategy for the assessment of the "druggability" of STAT3 dimers and allow investigations of reported activating and inhibiting STAT3 mutants at the atomistic level of detail. In the present study, we focused on the effects exerted by reported STAT3 mutations on the protein structure, dynamics, DNA-binding, and dimerization, thus linking structure, dynamics, energetics, and the biological function. By employing atomistic molecular dynamics and umbrella-sampling simulations to a series of human STAT3 dimers, which comprised wild-type protein and four mutations, we explained the modulation of STAT3 activity by these mutations. Counter-intuitively, our results show that the D570K inhibitory mutation exerts its effect by enhancing rather than weakening STAT3-DNA interactions, which interfere with the DNA release by the protein dimer and thus inhibit STAT3 function as a transcription factor. We mapped the binding site and characterized the binding mode of a clinical candidate napabucasin/BBI-608 at STAT3, which resembles the effect of a D570K mutation. Our results contribute to understanding the activation/inhibition mechanism of STAT3, to explain the molecular mechanism of STAT3 inhibition by BBI-608. Alongside the characterization of the BBI-608 binding mode, we also discovered a novel binding site amenable to bind small-molecule ligands, which may pave the way to design novel STAT3 inhibitors and to suggest new strategies for pharmacological interventions to combat cancers associated with poor prognosis.

Introduction

Signal transducer activator of transcription 3 (STAT3) protein has emerged as a prominent target in tumor progression due to its pivotal role in cell signaling.[1] The activation of the STAT3 protein has been also related to drug resistance,[2] to the expression of other anti-apoptotic proteins,[3] and to the inflammatory processes in tumor development, among others.[4−6] Despite its importance in cancer progression, the pharmacological drugging of STAT3 is still a challenge that demands clarification. Its tendency to aggregate is a major hurdle that prevents the determination of the structure in both monomeric and dimeric forms as well as bound to small-molecule inhibitors.[7−9] Although many strategies have been described in the literature to inhibit STAT3, only a few inhibitors are still going through clinical trials (e.g., TTI-101 [ClinicalTrials.gov Identifier: NCT03195699] or napabucasin (BBI-608)[10−12] [ClinicalTrials.gov Identifier: NCT03647839]), and STAT3 has become one of the most challenging cancer-related proteins to target by small molecules.

Gaining insights into the atomistic-level characteristics of STAT3 would permit the identification of novel strategies to interact with this protein by small molecules and target oncogenic pathways indirectly. The human STAT3 monomer is composed of six highly specialized domains (i.e., N-terminal, coiled-coil domain, DNA-binding domain, linker domain, SH2 domain, and C-terminal transactivation domain). Particularly, the DNA-binding domain (residues: 320–494) is responsible for the DNA-binding when STAT3 is in the dimeric form. An α-helix linker domain joins the latter with the SH2 domain, which is essential for the binding of STAT3 to phosphorylated receptors and for its dimerization (residues: 493–583). This process is mainly facilitated by the SH2 domain, as each monomer interacts after the phosphorylation of a specific tyrosine residue (Y705) located in the transactivation domain.[13,14] Therefore, in an attempt to avoid the dimerization, the SH2 domain has traditionally been the main target for drug design, mostly accompanied by computational studies relying on molecular docking calculations or similar structure-based approaches,[15−23] despite no crystallographic data being available up to date to support them. These ligands attempt to compete with phosphorylated p-Y705 at the site known to recognize phosphorylated residues[24,25] with limited success. An alternative, represented by OPB-31121[26] and OPB-51602,[4] is to target an allosteric site at the SH2 domain: these compounds bind to a pocket different from the one binding p-Y705.[4,26] Furthermore, STAT3 can undergo other post-translational modifications besides Y705 phosphorylation, such as S727 phosphorylation,[13,27] and it has been experimentally demonstrated that unphosphorylated STAT3 can dimerize and bind to DNA. This provides an alternative strategy to targeting the SH2 domain directly. Due to the challenges coupled with the effective targeting of STAT3 SH2 domains, other STAT3 domains have thus been explored in the development of potent and selective STAT3 inhibitors. Recently, several inhibitors have been identified to bind to the DNA-binding domain (DBD), which has been initially overlooked due to its potential specificity problems. Experimental evidence shows that its binding by small molecules permits the dimer formation, whereas disrupts the DNA–protein binding.[28−30] The exploration through alternatives to the SH2 domain might be the key to unlock the STAT3 “druggability” problem. The conservation of the three-dimensional structure of the whole protein is crucial for the activity of STAT3. Experimental data have demonstrated that point mutations in the linker domain suggest contacts with both the DNA-binding and SH2 domains, which could cause structural changes that severely affect STAT3 activity.[31] Alanine scanning demonstrated that the modification of interdomain hydrogen bonds can produce a significant decrease (i.e., K551A, W546A) or increase (D570K) in the STAT3–DNA-binding compared to that in the wild-type protein.[31] Understanding the effect of point mutations on STAT3 activity at the atomistic level could provide significant information about novel binding sites, unveiling new ways to target STAT3 by small-molecule ligands.

The lack of a ligand-bound crystallized structure hinders the path for an effective structure-based ligand design for both computational and medicinal chemists. Only a few moderate-resolution crystal structures of mouse STAT3 (PDB id: 1BG1,[9]3CWG,[7] and 4E68(8)) are available, and none of them present a bound ligand. A small number of STAT3 inhibitors are under clinical trials, but in most cases their binding sites at STAT3 and hence their binding modes remain unknown. Among them, napabucasin (BBI-608) is a first-in-class cancer stemness inhibitor that targets STAT3,[12] which is being tested (phase 3) as a treatment in advanced colorectal cancer.[11] Ji and co-workers reported that napabucasin binds in a pocket between the linker and the DNA-binding domain in a STAT3 crystal structure,[32] but the crystal structure has not been disclosed.

In this study, we addressed the druggability of the STAT3 dimer at the atomistic level of detail by performing equilibrium all-atom molecular dynamics (MD) and umbrella-sampling (US) simulations to determine the behavior of the wild-type and mutated STAT3. Monitoring the time evolution of the different domains at a time scale of hundreds of nanosceconds revealed their differential behaviors in terms of molecular flexibility and their ability to bind DNA when point mutations were introduced. The examination of mutation-induced structural changes directed us to new sites amenable to inhibition by small-molecule ligands. Furthermore, the construction of a homology model of ligand-bound STAT3 along with MD and US simulations helped us understand the behavior of a ligand that interacts with a newly identified “druggable” region of STAT3. Since this study focuses on the DBD and its interaction with DNA, we took special interest in Husby’s previous work, which investigated the interaction between STAT3 and DNA.[33]

Results and Discussion

Mutations Directly Affect the Sampled Configurations of the STAT3 Dimer

In the study by Mertens and co-workers, several mutations were indicated as crucial to control the DNA retention time within its respective binding cleft at STAT3.[31] These mutations occurred either in the DNA-binding domain or in the interdomain region. Prior to the US simulations to evaluate the DNA-binding in the mutated systems, a series of 50 ns MD simulations are performed to equilibrate the system and identify different behaviors of the studied systems with regard to the WT. Systems have been studied for the amount of time that Husby’s work claims to be necessary to achieve an energetically conserved and stable simulation.[33] Our results agree with those of Husby, as we could observe how one of the protein monomers is more flexible than the other and that root-mean-square deviation (RMSD) variations are mainly capitalized by the loops of SH2 and the coiled-coil domains.

One of the most significant configurational changes occurred within the D570K mutant, as the DNA double helix shifted downward (Figure ). This was most likely caused by the electrostatic effects at the residue located in the interface between the linker and the DNA-binding domain. The modification of the side-chain charge from negative (D) to positive (K) increased favorable protein interaction with the negatively charged nucleic backbone. This tightened the DNA-binding, resulting in a higher average DNA RMSD when compared to that of the crystal structure. It strongly indicates that the end-point configurations of the protein–DNA complexes play a significant role in their binding free energy, since the protein–nucleotide interactions change significantly between different mutants.

Figure 2

MD simulations show conformational changes between WT (blue) and D570K (red) STAT3 dimer systems. In the D570K mutant, one of the monomers is shifted (panel A to panel B), changing the conformational landscape of the dimer. Panel C shows how the position of the DNA duplex is shifted downward in the D570K mutant compared to the WT.

We assessed whether the conformational changes induced by the D570K mutation were observed in other mutations. Figure shows root-mean-square deviation (RMSD) plots of all STAT3 considered in this study as well as the WT protein. Except D570K, there were no large differences in protein RMSDs between the mutated STAT3 dimers and WT. This gap between D570K and other mutants is likely to arise from the combination of electrostatic and steric effects (all other mutations replaced large and polar residues with the smaller and apolar alanine), which affect the intrinsic dynamics of the coiled-coil (CC) domain. Hence, the dynamics of the CC domain might “tune” DNA–STAT3 interactions by allowing adjacent SH2 and DBD domains to improve their structural “fit” to the DNA.

Figure 3

Protein (A) and DNA (B) RMSDs after 50 ns of MD simulation. The D570K (blue) mutation shows a higher RMSD in both protein and DNA counterparts compared to the other mutations and WT-STAT3.

To follow up on the effects of the mutations that control DNA retention at the STAT3 on the structure, dynamics, and energetics of STAT3–DNA complexes, we carried out a set of umbrella-sampling (US) simulations, where DNA has been pulled from the STAT3 dimer. This process was studied using a series of US windows and simulated for 10 ns each.

The potential of mean force (PMF) calculated via the weighted histogram analysis method (WHAM) was consistent with the results reported by Mertens and co-workers[31] in most cases. Experimental results showed a drop in the DNA-binding for the tested interdomain mutations (EE434/435AA, W546A, and K551A) through time and an extraordinarily high retention time for D570K, with 100% DNA-binding even after 2 h.[31] The WT STAT3 had a higher energy barrier to reach its unbound state in comparison to the W546A mutant (Table ). The mutant shows a lower retention time, which indicates that interdomain interactions between mutated residues and E434 are crucial for STAT3–DNA-binding. Therefore, disrupting these interactions could represent an attractive strategy to target STAT3 by small molecules. K551A shows a profile similar to that of D570K, but presents a PMF value lower than that of the mutation and a few kcal/mol more than that of the WT (Figure and Table ).

Table 1

Free Energy Change Calculated via US for all US Calculations

run	ΔGUS (kcal/mol)
WT	–42.6 ± 2
EE43435AA	–62.7 ± 2
W546A	–35.7 ± 2
K551A	–58.7 ± 1
D570K	–58.5 ± 1
BBI-608	–127.0 ± 1

Figure 4

Potential of mean force (PMF) of dissociation of DNA from the STAT3 dimer. The W546A (purple) mutant shows a lower PMF than the WT (violet), whereas the PMF of K551A (brown) mutant is in the same range as that of the WT. Both EE434435AA (green) and D570K (marine green) would have displayed similar results, but the interaction between the DNA duplex and DBD of one STAT3 monomer resulted in higher PMF values than expected in later sampling windows. This indicates that DNA–protein interactions are more favorable in these mutants relative to the WT. The BBI-608 binding (blue) showed a trend similar to that of EE434435AA and D570K mutations, but much more pronounced.

The equilibrium MD simulations as well as PMF showed that the binding affinity of D570K to DNA is more favorable than that of any other mutant and more than that of WT STAT3 (Figure ). This indicates that this mutation promotes a tight binding between STAT3 and DNA, with a higher energy gap for DNA release upon pulling. The PMF curve showed that DNA-pulling from D570K required a higher energy gap to release the DNA from the STAT3 dimer. The experimental retention time correlated with the simulations when compared to WT STAT3. The data showed that DNA-binding to D570K was persistent through time, did not drive transcription, and resisted dephosphorylation, and thus prevented STAT3 from exerting its function.[12]

Collectively, these results indicate that D570K promotes a tight DNA-binding, so much that the bound duplex stays “locked” between the dimers, which effectively inhibits STAT3 by preventing it from releasing DNA and exerting its function as a transcription factor.

The major discrepancy between our simulations and experimental data was observed for the EE434/435AA double mutant. In the simulations, the mutant showed a higher PMF value than the WT, which indicated that its DNA-binding affinity should be higher, whereas experimental data showed that its behavior resembled that of the W546A mutant, which has shown a considerable drop of DNA-binding through time. Analysis of the final US windows indicated that the middle of the DNA duplex interacted favorably with the DNA-binding domain of one STAT3 monomer, but not another. Therefore, the US curve of the EE434435AA mutant displayed higher values arising from these interactions (DNA–STAT3 monomer) rather than from favorable interactions of the DNA–STAT3 dimer, as was the case with the D570K mutant.

Mertens’ results[31] are based on 60 to 100 min experiments that evaluate the percentage of DNA-binding of STAT3. This process would include several STAT3 molecules that would most likely go through several mechanistic cycles. Therefore, it is plausible that we did not sample the conformations that are contributing to these results.

Arginine R414 Acts as a “Gatekeeper” for DNA-Binding

All simulations indicated significant conformational changes of the arginine R414, which were required to release DNA (or to allow the DNA-binding to the STAT3 dimer). R414 has been shown previously as one of the main residues for DNA identification and binding,[33] but was not disclosed as a gatekeeper residue. It is at a close distance from the DNA, and its initial position did not allow the DNA to exit from the dimer. By acting as a gatekeeper of DNA-binding, R414 exerted a key role in controlling the opening and closing of the STAT3 dimer as well as tuned the dynamics of STAT3 monomers by modulating intradomain DBD–SH2, CC–SH2, DBD–LD, and CC–LD interactions.(Figure )

Figure 5

Conformational changes of arginine R414 observed during our simulations. Evolution of the DNA duplex and R414 position over the simulation time is depicted by colors, from red to white. Along the DNA-pulling pathway from the STAT3 dimer, the R414 side-chain rotates, allowing the dissociation to occur.

Inhibition of STAT3 by Napabucasin (BBI-608)

The results of simulations of apoSTAT3 (WT and mutants) highlighted a set of interdomain residues, explaining their effect on STAT3 behavior and function at the atomistic level of detail. These observations may pave the way to novel strategies for STAT3 inhibition using small-molecule ligands. Traditionally, the structure-guided ligand design for STAT3 inhibition focused on targeting the SH2 domain, which intended to prevent STAT3 dimerization. Although several inhibitors have been described to bind to the SH2 domain,[18,34−38] there is no crystal structure available to confirm it, and most studies focusing on the SH2 domain binders relied mainly on semi-rigid molecular-docking calculations, which may raise the question as to whether the SH2 domain is the binding site for their studied molecules. A few inhibitors such as OPB-31121 and OPB-51602 proved their binding to an allosteric site in the SH2 domain via mutagenesis analysis,[26] giving some insight into an alternative pocket to study the SH2-mediated inhibition. Additionally, there are reports of a few inhibitors that target DBD. Their binding modes have been virtually assessed, but again, these studies rely heavily on molecular docking.[28−30]

Recently, Ji and co-workers reported that napabucasin (BBI-608), which is one of the few STAT3 inhibitors in advanced clinical trials (phase 3), binds to a small pocket between the linker and the DNA-binding domain in a STAT3 crystal structure.[32] Since the crystal structure has not been released to the public domain, we have assessed the druggability of this segment of STAT3, identified the putative pocket, and subsequently built the model of BBI-608 bound to STAT3 using the molecular-docking approach, and then validated the obtained binding mode by the atomistic MD and US simulations.

Molecular docking was performed by MOE for each STAT3 monomer separately, trying to generate the most plausible conformation relying on the limited data available. Both blind (whole monomer) and targeted (residues of the identified pocket) docking calculations resulted in a set of conformations with favorable energy scores and a highly populated cluster located within the DBD site pocket, in close contact with residues H332, P333, R335, K573, and D570 (Figure ). Two conformations, matching the published data, were found; both were assessed and validated.

Figure 6

Binding pose of BBI-608, deduced from the data reported by Ji and co-workers.[32] Since the crystal structure has not been released, we have modeled the most plausible binding mode by molecular docking. The side-chain of the mutated K570 residue is displayed and colored green; it is overlapping with the plausible location of BBI-608.

To validate the binding mode of BBI-608, MD simulations of the WT STAT3–DNA–BBI-608 complex were performed for 100 ns in triplicate. Subsequently, the ligand affinity has been calculated. The ligand docked on either of the SAT3 monomers remained bound through the whole simulation. Interaction energies calculated by MMPBSA analysis (g_mmpbsa[39] module) resulted in −18.1 ± 2.6 kcal/mol, showing a favorable binding.

US simulations, performed using the same protocol as for STAT3–DNA complexes, started the pull from the most populated cluster. The calculated binding affinity has been overestimated (−160 kcal/mol); nevertheless, it showed a very tight binding. Compared with the results obtained for protein–DNA complexes described in the previous section, it implies that the presence of BBI-608 enhances DNA-binding with an effect similar to that of the D570K mutation. As such, BBI-608 inhibits the function of STAT3 in a manner similar to the D570K mutation, which does not drive transcription and resists phosphorylation.[12] Since D570 has been annotated by Ji et al.[32] as the BBI-608-binding-site residue, we concluded that BBI-608 binding to WT STAT3 generated a DNA–protein interaction pattern and retention time similar to that of the D570K mutation.

BBI-608 Binding Promotes STAT3 Dimerization

We performed MD simulations of the BBI-608-bound WT STAT3 dimer without DNA to study the influence of the ligand on the protein behavior in the absence of DNA. BBI-608 was bound to each STAT3 monomer, and four 100 ns replicas showed a variation of results. In only one out of the four replicas, both BBI-608 molecules remained bound in their pockets through the whole trajectory, whereas DBD domains of both STAT3 monomers moved closer to each other, reaching the point of forming interdomain hydrogen bonds. In one simulation, both ligand molecules dissociated from the pockets, which caused an opening of the STAT3 dimer, and in the other two, only one of the ligands left the binding cavity at the 35 and 70 ns marks.

To assess ligand-induced conformational changes within the DNA entry through both DBD domains, distances between some of the residues involved in H-bonding (e.g., Q344-G342 and T412-Q344) were measured and analyzed in all four replicas and compared to the simulation of the WT STAT3 dimer without DNA and/or BBI-608 bound (Figure ).

Figure 7

(A) Interatomic distances between the Cα of residues Q344 and G432 and residues T412 and Q344, calculated along a 50 ns MD simulation of ligand–STAT3 complexes. Replica 1 (blue) consists of the dissociated system and both Replicas 2 and 3 (orange and green) keep their ligand bound through the whole simulation. (B) Close-up and (C) panoramic views as the STAT3 dimer closes once the BBI-608 molecule is bound. (D) Protein–ligand energy interaction for both BBI-608 molecules interacting with each monomer along a 50 ns MD simulation (three replicas: blue, red, and green). The STAT3 dimer bound to DNA (black) is shown as the reference.

The STAT3 dimer closed further down in the presence of BBI-608. This was particularly pronounced in one of the replicas, in which the distance between monomers reduced to <5 Å. These results indicate that BBI-608 binding to apoSTAT3 is likely to trigger conformation changes that would prevent DNA from binding. In the replica simulation, where both ligands dissociated from STAT3, the distance between monomers increased upon ligand dissociation, as both ligands exit via the gap formed between DBD domains of STAT3 monomers.

The protein–ligand interaction energy was calculated, and a correlation between ligand dissociation, dimer separation, and poor ligand affinity could be observed.

The simulations also indicate that the ligand binding to one of the STAT3 monomers is more favorable than it binding to another one. Although interaction energies are favorable for both monomers, one shows an interaction energy value twice as favorable as for another monomer. Although the allosteric effects within STAT3 were beyond the scope of this study, these results strongly suggest that such effects may occur in STAT3 dimers and contribute to the modulation of STAT3 by inhibitors. Another theory could be that the model is not as optimal as expected. Although the used docking poses have an extraordinary resemblance to Ji’s data,[32] only one is shown in the patent. Since this dimer is asymmetrical, DNA interactions with both cavities will be different and could imply a different conformation and/or location for the ligand that does not correspond to the model one. Furthermore, the MD results strongly indicate that one monomer is much more mobile than the other, and this is likely to affect protein–ligand affinity fluctuations.

Identification of a Novel Druggable Binding Site

Experimental results[12] combined with the simulations strongly indicated that targeting the interface between STAT3 monomers may trigger a response similar to the inhibition by ligands binding to the DBD domain, and therefore must be explored in structure-based ligand design efforts. With most STAT3 ligands being designed for the SH2 domain and just a few for the DBD domain, the identification of new druggable pockets for STAT3 inhibition is of great interest.

As such, we scanned STAT3 for the presence of potential binding sites with two pocket detection tools: fpocket[40] and MOE’s Site Finder. Upon selection of the dimer model and different clusters from its MD simulation trajectories, we confirmed the binding site identified for BBI-608,[32] which has been identified by both tools as their top-ranked site. Interestingly, this site was detected by both fpocket and SiteFinder for all analyzed structures (Figure ). In addition, a new pocket within the DNA-binding domain (in close contact with E434 and E435) was identified. The main difference between results obtained by both tools is that fpocket was more prone to detect SH2 domain sites as pockets (Figure , B2-3), whereas SiteFinder identified a novel druggable pocket close to R414 (Figure , A1). Ligand binding to the R414 pocket could result in a possible DNA-release/binding impediment.

Figure 8

Druggability of the STAT3 dimer. Sitefinder (A) and fpocket (B) were used to identify new potential pockets for structure-based drug design. In both cases the BBI-608 DBD site was identified along with novel DBD pockets.

Conclusions

The use of computational approaches based on atomistic molecular dynamics simulations enabled us to understand the effects of specific STAT3 mutations, which were described in the literature, and to explain their modulation of the STAT3 activity. Consistent with the findings of Mertens and co-workers,[32] D570K mutation exerted its effect by enhancing interactions between STAT3 and DNA, which interfered with DNA release by the STAT3 dimer and thus inhibited the protein’s function by not driving transcription and resisting dephosphorylation. Subsequently, recent identification of a plausible binding site for the small-molecule STAT3 inhibitor nababucasin (BBI-608) helped us to deconvolute its inhibition mechanism, which resembled the effect exerted by the D570K mutation. The identification of the putative binding site for BBI-608 at the DNA-binding domain may contribute to the development of novel potent and selective STAT3 inhibitors. The similarity between the model and the available patent data[32] raise the question as to why this pocket has been overlooked before. Mutated interdomain residues E435, W546, and K551 unveil a poor binding to DNA, leading to yet another way of targeting STAT3 by disrupting theses residues, pointing toward possible novel allosteric binding sites. Structure-based ligand designs targeting these novel pockets, coupled with novel methodologies, such as employing the recently developed FragLites,[41] are likely to expand a set of chemotypes active toward STAT3 and contribute to the development of novel inhibitors of this important yet very challenging drug target.

Methods

Molecular Modeling of the Human STAT3 Dimer

Initial models of dimeric human STAT3 (wild type and mutants) in a complex with DNA were created using the crystal structure of unphosphorylated mouse STAT3B (PDB code: 4E68), which spans residues 136–716. The N-terminal domain has been excluded from the structures subjected to the simulations. Loops spanning the residues 184–194 and 688–702 were modeled using the MODELLER interface[42,43] available in UCSF Chimera.[44] The DNA double-strand bound to the model was designed based on 4E68, with the 5′–3′ strand sequence of TGCATTTCCCGTAATC. The final model was subjected to 20 000 cycles of the steepest descent energy minimization.

The STAT3 mutations (Figure ), which were selected following the study by Mertens and co-workers,[31] were introduced in UCSF Chimera by swapping side-chains to the target residues and adjusting new conformations using Dunbrack’s rotamer library integrated within UCSF Chimera.

Figure 1

View of the mutated residues that form interdomain interactions and affect STAT3 (PDB code 1BG1) behavior, particularly DNA-binding. SH2 domain (cyan), linker domain (orange), DNA-binding domain (red), and DNA (blue) are highlighted. Side-chains of the target residues are displayed as gray sticks.

Modeling of Ligand-Bound STAT3

The crystal structure of ligand-bound STAT3 is not yet available; therefore, we built the most similar model possible with the accessible data. Starting from the dimer model of wild-type STAT3 described in the previous section, molecular-docking calculations were performed with MOE,[45] using napabucasin (BBI-608) as the ligand. To ensure the scoring function accuracy with the target and the best possible fit, both blind (full dimer, residues 136–716) and targeted docking (pocket described) were performed. Two hundred different conformations of the ligand were scored for each run using Triangle Matcher and London ΔG for the first scoring function. Thereafter, the top 100 conformations were rescored using Induced Fit and the generalized Born volume integral/weighted surface area score.[45] From the final poses obtained, one for each monomer was selected based on the score, interactions, and consistency with the experimental structure shown in Ji’s work,[32] which shows BB1-608 bound to a pocket located within the DNA–STAT3 interface.

Molecular Dynamics and Umbrella-Sampling Simulations

To study the effects of the dynamics of mutations in the STAT3 dimer, each of the listed mutations was created in UCSF Chimera.[44] Structural hydrogens were added and the following protein parametrization was performed using the Gromacs 2016.03[46] suite with AMBERFF99SB-ILDN[47] force field. Before pulling the DNA from the complex, the systems were relaxed with a short equilibrium MD production run. Hence, a 1 nm cubic box was centered on the structure and the system is solvated with TIP3P waters. Sodium and chloride ions were added to a concentration of 0.1 M, resulting in systems with more than seven hundred thousand atoms. Bonds were constrained using the LINCS[48] algorithm. The electrostatic interactions were calculated using the particle-mesh Ewald method, with a nonbonded cut-off set at 0.1 nm. All structures were minimized via the steepest descent algorithm for 20 000 steps of 0.02 nm, and minimizations were stopped when the maximum force fell below 1000 kJ/(mol nm) using the Verlet cut-off scheme.[49] After minimization, temperature equilibration was performed for 100 ps with a time step of 2 fs with position restraints applied to the backbone using an NVT. The temperature coupling was set between the protein and the nonprotein entities by using a Berendsen thermostat,[50] with a time constant of 0.1 ps, and the temperature set to reach 300 K with the pressure coupling off. Sequentially, a pressure NPT equilibration was performed, followed by 100 ps of NVT equilibration, 100 ps of NPT equilibration, and a production run of 100 ns. The temperature was set constant at 300 K by using a modified Berendsen thermostat (τ = 0.1 ps).[50] The pressure was kept constant at 1 bar by Parinello–Rahman isotropic coupling (τ = 2.0 ps) to a pressure bath.[51]

For the umbrella-sampling simulation, a 50 ns pre-umbrella equilibration was carried out, with the complex rotating its principal axis to align with the z-axis of the simulation box. A pull sampling was performed using a constant force approach (k = 1000 kJ/(mol nm), with a rate of 0.01 nm) between the centers of masses of the SH2 domain and the DNA double helix, along the described path shown in Figure . From each corresponding pull simulation, a series of conformations have been selected to sample the process of entering–exiting the DNA-binding site. Each of the 25 selected umbrella windows has been through a 1 ns NPT equilibration run, followed by a 5 ns NPT distance-restrained production run, totalizing per system 135 ns of simulation time, using the previously described protocol and parameters. Afterward, the potential of mean force (PMF) curve of the studied scenario has been calculated with the weighted histogram analysis method (WHAM) tool from Gromacs,[52] and associated errors were calculated using both the convergence criteria and the implemented bootstrap method in Gromacs WHAM. All calculations for the analysis were made using the Gromacs tools.

Figure 9

Description of the STAT3 umbrella-sampling (US) simulations performed in this study. After 100 ns of equilibrium MD production run, DNA is pulled out of the dimer and a series of windows are selected that will be simulated for 10 ns to compute the potential of mean force (PMF) of DNA.