Targeting SHP2 Cryptic Allosteric Sites for Effective Cancer Therapy

SHP2, a pivotal component downstream of both receptor and non-receptor tyrosine kinases, has been underscored in the progression of various human cancers and neurodevelopmental disorders. Allosteric inhibitors have been proposed to regulate its autoinhibition. However, oncogenic mutations, such as E76K, convert SHP2 into its open state, wherein the catalytic cleft becomes fully exposed to its ligands. This study elucidates the dynamic properties of SHP2 structures across different states, with a focus on the effects of oncogenic mutation on two known binding sites of allosteric inhibitors. Through extensive modeling and simulations, we further identified an alternative allosteric binding pocket in solution structures. Additional analysis provides insights into the dynamics and stability of the potential site. In addition, multi-tier screening was deployed to identify potential binders targeting the potential site. Our efforts to identify a new allosteric site contribute to community-wide initiatives developing therapies using multiple allosteric inhibitors to target distinct pockets on SHP2, in the hope of potentially inhibiting or slowing tumor growth associated with SHP2.


Introduction
The SHP2, a non-receptor protein tyrosine phosphatase encoded by PTPN11, is integral to the developmental processes mediated by growth factors, cytokines, and adhesion receptors.SHP2 is ubiquitously expressed and is essential for the sustained activation of the Ras-MAP kinase pathway, which plays a crucial role in cell proliferation, differentiation, and survival.Additionally, SHP2 modulates several other critical signaling pathways, including NF-κB, JAK-STAT, PI3K-AKT, PD-1, and various immune checkpoints, which are involved in immune regulation and cell growth [1][2][3][4].
Activating mutations in SHP2 are implicated in developmental disorders such as Noonan Syndrome [5] and are frequently observed in juvenile myelomonocytic leukemia (35%), and, to a lesser extent, in acute myeloid leukemia (5%) [6].These mutations are also present at lower frequencies in other hematologic and solid tumors [4] and have been shown to induce leukemia in murine models [7].Conversely, inhibiting SHP2 demonstrates antitumor activity across various cancer models [8,9].
Int. J. Mol.Sci.2024, 25,6201 2 of 15 C-SH2 domain, represents a conformational change of unusual magnitude that has garnered significant research attention [11][12][13][14].Nuclear magnetic resonance (NMR) spectroscopy and molecular dynamics simulations have been instrumental in elucidating the key transition states and pathways underlying SHP2's allosteric regulation [12,[15][16][17].Recent investigations, employing molecular dynamics simulations and small-angle X-ray scattering (SAXS), have revealed unexpected flexibility within the open state of the E76K mutant, a model system for studying SHP2 activation [16,[18][19][20].Notably, molecular dynamics simulations show that the E76K solution structure deviates significantly from the crystal structure [18].Furthermore, Anselmi and Hub reported a heterogeneous atomistic ensemble of the E76K mutant in solution, in excellent agreement with SAXS data [20].These observations raise questions about the factors contributing to this flexibility, including potential intrinsic domain instability or the influence of crystal packing effects.
The development of allosteric inhibitors that stabilize the closed state of SHP2 has shown promise in clinical trials for advanced or metastatic solid tumors [21][22][23][24][25][26][27][28][29][30].These "molecular glue" inhibitors, like SHP099, target an allosteric site formed between the C-SH2 and PTP domains, demonstrating high-affinity binding and antiproliferative effects in certain cancer cell lines [8].However, activating mutations can destabilize the closed state, leading to reduced inhibitor efficacy and drug resistance [31].While a dual inhibition strategy targeting an additional allosteric site at the N-SH2/PTP interface has been proposed [32], the dynamic nature of SHP2 poses challenges for inhibitor binding, highlighting the need for more robust allosteric inhibitors.
In light of these challenges, this study aims to elucidate the origin of the observed flexibility within the open state of SHP2, particularly in the context of the E76K mutant.We will also investigate the impact of activating mutations on the efficacy of existing allosteric inhibitors.These insights will contribute to the development of more effective therapeutic strategies targeting SHP2, including the identification of novel allosteric sites less susceptible to the effects of activating mutations.Recent investigations, employing molecular dynamics simulations and small-angle X-ray scattering (SAXS), have revealed unexpected flexibility within the open state of the E76K mutant, a model system for studying SHP2 activation [16,[18][19][20].Notably, molecular dynamics simulations show that the E76K solution structure deviates significantly from the crystal structure [18].Furthermore, Anselmi and Hub reported a heterogeneous atomistic ensemble of the E76K mutant in solution, in excellent agreement with SAXS data [20].These observations raise questions about the factors contributing to this flexibility, including potential intrinsic domain instability or the influence of crystal packing effects.
The development of allosteric inhibitors that stabilize the closed state of SHP2 has shown promise in clinical trials for advanced or metastatic solid tumors [21][22][23][24][25][26][27][28][29][30].These "molecular glue" inhibitors, like SHP099, target an allosteric site formed between the C-SH2 and PTP domains, demonstrating high-affinity binding and antiproliferative effects in certain cancer cell lines [8].However, activating mutations can destabilize the closed state, leading to reduced inhibitor efficacy and drug resistance [31].While a dual inhibition strategy targeting an additional allosteric site at the N-SH2/PTP interface has been proposed [32], the dynamic nature of SHP2 poses challenges for inhibitor binding, highlighting the need for more robust allosteric inhibitors.
In light of these challenges, this study aims to elucidate the origin of the observed flexibility within the open state of SHP2, particularly in the context of the E76K mutant.We will also investigate the impact of activating mutations on the efficacy of existing allosteric inhibitors.These insights will contribute to the development of more effective therapeutic strategies targeting SHP2, including the identification of novel allosteric sites less susceptible to the effects of activating mutations.

Results and Discussion
In this study, we employed molecular dynamics (MD) simulations to investigate the solution structures and dynamic stability of SHP2 in both its wild-type (closed state, wtSHP2) and E76K mutant (open state, mtSHP2) forms.Simulations were conducted using the FF19SB and GAFF2 force field in the TIP3P water and 298 K. Three production trajectories of 500 ns each were collected on both monomeric (chain A and chain B) and dimeric (chain A + chain B) configurations.To assess the impact of ligand binding, we also simulated systems with the allosteric inhibitor SHP099 (chain A), also in three replicas of 500 ns each.Initial structures for all simulations were derived from available SHP2 crystal structures.A total of approximately 28.6 µs of MD trajectories were collected in the explicit water model (see Supplementary Table S1 for details).
Our analysis first focused on the dynamic stability of individual domains (N-SH2, C-SH2, and PTP) to quantify the increased flexibility observed in the E76K open state.As detailed in the Supplementary Information (Figure S1), each domain within wtSHP2 exhibited low RMSD values, indicating remarkable stability.This is likely due to the compact, tightly packed arrangement of the protein in the closed state.In contrast, while the N-SH2 and PTP domains of mtSHP2 maintained RMSD values comparable to wtSHP2, the C-SH2 domain displayed significantly higher fluctuations, suggesting heightened flexibility.This observation aligns with the lower certainty of the C-SH2 domain, which reflects its relatively looser packing environment in the crystal lattice, both intermolecularly and intramolecularly.
We further explored the role of crystal packing environment, a notable factor influencing structural variations in MD simulations.While four monomers are organized linearly in a relatively straightforward packing pattern within the wtSHP2 unit cell, eight monomers arrayed in a more complex packing formation within the mtSHP2 unit cell (Supplementary Information, Figure S2).Obviously, the wtSHP2 unit cell closely mirrors our solvated monomer/dimer MD simulation environment.In contrast, the mtSHP2 asymmetric unit exhibits a more complex packing arrangement within the unit cell, leading to solvated monomer/dimer MD simulations that do not fully capture this intricate packing.Thus, the high variability in C-SH2 domain dynamics in mtSHP2 simulations may be partly due to the absence of crystal packing constraints.Conversely, the lower C-SH2 variability in wtSHP2 simulations could be attributed to solvation effects that better mimic the wtSHP2 crystal environment.
Given that the observed variations in individual domain dynamics are primarily attributed to differences in crystal packing, the success of allosteric modulators in shifting the open state structure towards the closed state is reinforced.The primary role of these modulators is to act as "molecular glues", stabilizing inter-domain interactions.

Impact of the Oncogenic Mutation E76K on Allosteric Inhibition Sites
Our study extends to the functional disruption of the E76K mutation on SHP2's allosteric inhibition.Traditionally, allosteric SHP2 inhibitors have focused on the allo-site-1 pocket [8], which is formed by residues from both the C-SH2 and PTP domains, such as Arg111, Phe113, Glu250, Leu254, Gln257, Pro491, and Gln495 (Figure 2D,E).These inhibitors function by stabilizing SHP2 in its closed state [21][22][23][24][25][26][27][28][29][30].Through monitoring the CA distance between Arg111 (on C-SH2 domain) and Phe113 (on PTP domain) across various SHP2 states, we observed a uniform distance in all closed conformations, including wtSHP2-ub, wtSHP2-b, and mtSHP2-b (Figure 2C).However, the open state, influenced by the E76K mutation, exhibits a pronounced CA distance deviation, correlating with the C-SH2 domain's 120-degree rotation (Figure 2A,B) [11,12].This movement signifies the allo-site-1 pocket disruption upon SHP2 activation.In the closed state, SHP099 binds from the back side, while the front-side access is hindered by a short disordered loop.However, mutation E76K alters the binding pocket through a rotation in the C-SH2 domain, resulting in the displacement of key residues Arg111 and Phe113, which play a role in the creation of the binding pocket for SHP2 inhibitors.(C) represents the CA distances between Arg111 and Phe113 that are involved in SHP099 binding across all SHP2 configurations.(D,E) illustrate the superimposed structures from the top MD solution clusters of each SHP2 system, while the binding site residues form a consistent cavity in the wtSHP2 and mtSHP2-b states, a discernible shift in two residues from the C-SH2 domain away from the binding pocket is evident for the pocket disruption in the mtSHP2-ub conformation.The green sticks denote the crystallographic structure from the wtSHP2-b state.
A novel dual inhibition strategy has been proposed for more effective SHP2 inhibition, targeting a second pocket, allo-site-2, at the interface of the N-SH2 and PTP domains (Supplementary Figure S3) [32].Despite this, achieving potent inhibition remains challenging, indicating the potential for improved inhibitor design.We hypothesized that the reduced efficacy of inhibitors at pocket allo-site-2 could be due to its transient stability.The pocket is flanked by key residues Glu83 and Arg265, which may predispose the site to inhibitor dissociation.
To evaluate this hypothesis, we assessed the CA distance dynamics between Glu83 and Arg265, which are crucial to the structure of the allo-site-2 pocket.As shown in Figure 3A-D, these residues initially appear to be stable, maintaining the pocket's integrity.A novel dual inhibition strategy has been proposed for more effective SHP2 inhibition, targeting a second pocket, allo-site-2, at the interface of the N-SH2 and PTP domains (Supplementary Figure S3) [32].Despite this, achieving potent inhibition remains challenging, indicating the potential for improved inhibitor design.We hypothesized that the reduced efficacy of inhibitors at pocket allo-site-2 could be due to its transient stability.The pocket is flanked by key residues Glu83 and Arg265, which may predispose the site to inhibitor dissociation.
To evaluate this hypothesis, we assessed the CA distance dynamics between Glu83 and Arg265, which are crucial to the structure of the allo-site-2 pocket.As shown in Figure 3A-D, these residues initially appear to be stable, maintaining the pocket's integrity.However, the wild type complex (wtSHP2-b) gradually transitions into a disrupted state with a flat surface after ~300 ns, and the mutant complex (mtSHP2-b) transitions into a state with a salt bridge after ~100 ns, though occasional fluctuations back to the initial state are also visible.These conformational changes underscore the transient nature of allo-site-2's pocket stability.When these residues are close together, the ligand cannot bind.When these residues move apart, the pocket becomes a flat surface, both compromising the pocket's binding potential and resulting in the inhibitors showing less potency toward this site.
Int. J. Mol.Sci.2024, 25, 6201 5 of 15 However, the wild type complex (wtSHP2-b) gradually transitions into a disrupted state with a flat surface after ~300 ns, and the mutant complex (mtSHP2-b) transitions into a state with a salt bridge after ~100 ns, though occasional fluctuations back to the initial state are also visible.These conformational changes underscore the transient nature of allo-site-2's pocket stability.When these residues are close together, the ligand cannot bind.When these residues move apart, the pocket becomes a flat surface, both compromising the pocket's binding potential and resulting in the inhibitors showing less potency toward this site.An area adjacent to allo-site-2, comprised of the disordered loop spanning residues His85-Asp90 and Gln92-Val95, also appears to be significant in ligand binding (Figure 3E, highlighted in yellow).We observed this loop oscillating between two distinct states: An area adjacent to allo-site-2, comprised of the disordered loop spanning residues His85-Asp90 and Gln92-Val95, also appears to be significant in ligand binding (Figure 3E, highlighted in yellow).We observed this loop oscillating between two distinct states: an 'in' conformation that acts as a cap, potentially sealing off the pocket (wtSHP2-b), and an 'out' conformation that opens allo-site-2, facilitating inhibitor binding (mtSHP2-b).This loop's flexibility introduces a steric barrier that could deter the binding of a second inhibitor to allo-site-2, thereby influencing the overall inhibitor-enzyme interaction dynamics.
In summary, the structural dynamics of Glu83 and Arg265, along with the flexible loop adjacent to allo-site-2, play critical roles in the pocket's stability and accessibility.These factors must be considered in the design of more effective SHP2 inhibitors.
We postulate the presence of additional allosteric sites that could be harnessed to stabilize the closed state of SHP2, particularly to counteract the hyperactive E76K mutant's effects.Therefore, discovering a novel secondary binding pocket is a key objective of our research aimed at achieving potent and comprehensive inhibition of SHP2.
Integrating these findings, we delved into the dynamics of the allo-closite pocket, both with and without the SHP099 inhibitor, to confirm its stability in all closed states of SHP2 (Figure 4).Based on CA distances between contributing SH2 and PTP domain residues, our structural analysis revealed that the allo-closite pocket forms robustly in all closed SHP2 states, with constituent residues coalescing to establish a well-defined pocket.Conversely, in the mtSHP2 state unbound to SHP099, we detected dynamic convergence of the SH2 domains and divergent movement of the PTP domain residues.After a 250 ns simulation (Figure 4A), a notable CA distance divergence among the pocket-forming residues in mtSHP2 unbound to SHP099 indicated pocket destabilization (Figure 4B).Intriguingly, when bound to SHP099 (Figure 4C,D), the allo-closite pocket exhibited greater stabilitysurpassing even the apo wtSHP2-highlighting the potential exclusivity of our identified pocket in the closed state and its selective affinity, akin to that of SHP099.

Hit Compound Identification
Our newly identified pocket underwent virtual screening against a cancer-specific small molecule library from the NCI Open Database (https://cactus.nci.nih.gov/download/nci/index.html, accessed on 28 March 2023).This rigorous process involved a threetiered screening approach designed to pinpoint potential hit compounds.Following a consensus selection process after Tier I and Tier II screenings, the top 2000 compounds for both wtSHP2 and mtSHP2 receptors were determined (Supplementary Table S2), leading to the identification of 73 common top-hit compounds (Supplementary Table S3).Because of large amounts of polar and charged residues in the allo-closite and the pocket's shape, these top hits underwent further scrutiny by visual inspection of their interactions with pocket residues and their drug likeness, resulting in a selection of 18 promising, relatively planar candidates for detailed analysis (Table 1), hereafter referred to as theoretical binding affinity (t-BA) hits, as their selection was primarily based on their predicted t-BA values.

Hit Compound Identification
Our newly identified pocket underwent virtual screening against a cancer-specific small molecule library from the NCI Open Database (https://cactus.nci.nih.gov/download/nci/index.html,accessed on 28 March 2023).This rigorous process involved a three-tiered screening approach designed to pinpoint potential hit compounds.Following a consensus selection process after Tier I and Tier II screenings, the top 2000 compounds for both wtSHP2 and mtSHP2 receptors were determined (Supplementary Table S2), leading to the identification of 73 common top-hit compounds (Supplementary Table S3).Because of large amounts of polar and charged residues in the allo-closite and the pocket's shape, these top hits underwent further scrutiny by visual inspection of their interactions with pocket residues and their drug likeness, resulting in a selection of 18 promising, relatively planar candidates for detailed analysis (Table 1), hereafter referred to as theoretical binding affinity (t-BA) hits, as their selection was primarily based on their predicted t-BA values.
In pursuit of a broader scope of candidate compounds, we expanded our visual inspection beyond the top 2000 compounds to include those that did not reach the highest ranks.This comprehensive evaluation considered the drug likeness, the binding orientations, and interactions with critical residues across the three domains, recognizing that the molecule's alignment is pivotal for effectively stabilizing the domain interfaces.This exhaustive analysis yielded 21 additional potential compounds (Table 2).When contrasting these compounds, referred to as human-visualized interaction (h-VI) hits, with the initial set of t-BA hits (Table 1), no overlap was observed.Yet, a comparison with the 73 t-BA hits from Supplementary Table S3 unveiled five shared hit compounds, demonstrating some convergence between the selection methodologies.Next, we conducted an in-depth analysis of the protein-ligand interaction (PLI) profiles for the 39 compounds selected through both t-BA and h-VI criteria.This examination began with the top-hit compounds identified based on theoretical binding affinity (t-BA).The PLI profiles for wtSHP2-b (Supplementary Figure S4) revealed predominant interactions between the compounds and key pocket residues, including Asn10 and Glu15 from the N-SH2 domain, along with Asp106, Glu139, Pro144, Phe147, Cys174, Gln175, Glu176, and Leu177 from the C-SH2 domain.For mtSHP2-b (Supplementary Figure S5), the PLI profiles pinpointed critical interactions within the pocket, notably with His8, Glu15, Asn18, Arg23 from the N-SH2 domain; Asp106, Glu176 from the C-SH2 domain; and Asp241, Lys242 from the PTP domain.
In summary, nearly all 39 compounds meeting either t-BA or h-VI criteria form hydrogen bonds with multiple residues across different protein domains.Furthermore, many of these compounds contain aromatic rings, facilitating additional hydrophobic or cation-π interactions.

Potential Binding to Allo-Site-1
For comparison, we also conducted docking analysis of the 39 selected compounds against the allo-site-1 pocket (on the receptor structure without SHP099, as in Tier I).The analysis details and docking results are provided in Tables S4 and S5, respectively.Notably, none of the 39 compounds exhibited high affinity for allo-site-1, when compared with the binding affinity of −13 kcal/mol for SHP099 in the crystallographic pose [8].Upon visualization of the binding poses, only NSC-163300, NSC-637201, NSC-39918, and NSC-63675 were successfully inserted into the tunnel-like pocket.These compounds also displayed the lowest predicted binding affinities among the 39 compounds tested.This comparative analysis suggests that these compounds are less likely to compete effectively for binding to allo-site-1, supporting our goal of using the compounds in a dual allosteric inhibition scheme along with SHP099.

Stability Analysis via MD Simulations in Tier III Screening
Our Tier III screening involved MD simulations of complexes with selected hit compounds.We conducted a total of 78 MD simulations for wtSHP2and mtSHP2-ligand systems to assess the stability of these complexes.Each hit compound complex underwent a 100 ns production MD simulation, which allowed us to weed out weak binders within the allo-closite pocket.
To gauge their dynamic stability, we monitored the root-mean-square deviation (RMSD) of each hit compound.We established a threshold RMSD of 3 Å as the stability criterion; compounds exhibiting an RMSD less than this value were deemed stable and were selected for further scrutiny.Our visual analysis indicated that compounds with RMSD under 3 Å maintained significantly greater stability and consistent orientation within the pocket compared to those with RMSD values above 3 Å (Figure 5A,B).
Our findings identified 11 highly stable, 13 partially stable, and 15 unstable hit compounds for wtSHP2.Correspondingly, for mtSHP2, we discovered 11 highly stable, 12 partially stable, and 16 unstable compounds.Notably, the highly stable compounds demonstrated similar stability profiles in both wtSHP2 and mtSHP2, suggesting a favorable selectivity for the closed state of SHP2 (Figure 5C,D).The partially stable and unstable compounds, however, exhibited variability between the two receptor states.
Additionally, we extended the production molecular dynamic (MD) simulations for the 11 highly stable compounds up to 500 ns as the crystal structures and analyzed their RMSD fluctuations within allo-closite (Figure S8).Most compounds remain well-positioned in the pocket, with RMSD under 3 Å.These extended simulations further underscore the high stability of the identified compounds.
pounds were occasionally seen approaching SHP099 in the allo-site-1 pocket, induc significant rotation of the 1,4-dimethylpiperidin-4-amine moiety in SHP099.This o vation suggests the potential for crosstalk between the allo-closite and allo-site-1 poc possibly resulting in increased SHP2 stability through a dual cooperativity mechan Figure 5F depicts the 2D structures of select hit compounds that were consistently fo to target the closed state of the SHP2 protein.

Multi-Tier Molecular Docking
Given several closed states for SHP2 sharing similar structures with mutual RM less than 1 Å, we went ahead to choose different receptors at different stages in our d ing screening, including wtSHP2-ub (closed-conformation, PDB ID: 2SHP [10]) init and wtSHP2-b (closed-conformation, PDB ID: 5EHR [8]) and mtSHP2b (closed-co mation, PDB ID: 6CRG [11]) simultaneously at the second stage.
Prior to initiating our screening approach, we first utilized RDKit (rele 2023.03.1), an open-source cheminformatics software [34,35], to refine the 0.26 mi compounds from the NCI Open Database.Utilizing a customized Python script, we screened compounds by stripping metal-containing compounds, acknowledging thei tential off-target interactions and toxicity risks.We also set a molecular weight thres of 700 Daltons, emphasizing the importance of solubility and permeability in drug de Finally, our simulations revealed evidence of induced fit within the pocket upon ligand binding.The two critical gatekeeper residues, Arg23N-SH2 and Glu176C-SH2 (Figure 5E), initially positioned approximately 17 Å apart, were observed to form a robust salt bridge with a sub-3 Å distance upon full relaxation of the compounds at the end of the 100 ns MD, enhancing pocket stability.Additionally, the polar segments of certain compounds were occasionally seen approaching SHP099 in the allo-site-1 pocket, inducing a significant rotation of the 1,4-dimethylpiperidin-4-amine moiety in SHP099.This observation suggests the potential for crosstalk between the allo-closite and allo-site-1 pockets, possibly resulting in increased SHP2 stability through a dual cooperativity mechanism.Figure 5F depicts the 2D structures of select hit compounds that were consistently found to target the closed state of the SHP2 protein.
Prior to initiating our screening approach, we first utilized RDKit (released 2023.03.1), an open-source cheminformatics software [34,35], to refine the 0.26 million compounds from the NCI Open Database.Utilizing a customized Python script, we pre-screened compounds by stripping metal-containing compounds, acknowledging their potential off-target interactions and toxicity risks.We also set a molecular weight threshold of 700 Daltons, emphasizing the importance of solubility and permeability in drug design.Additionally, compounds with an absolute net charge of more than 1 a.u.are discarded to discourage artificially strong electrostatic interactions in docking and to ensure targeted binding efficacy and seamless membrane passage.
We begin with Tier I screening which includes docking of the pre-screened NCI compounds targeting the SHP2 basal state without the SHP099 inhibitor, i.e., wtSHP2-ub.The most dominant wtSHP2-ub MD solution conformation extracted from the k-mean clustering analysis was utilized as a receptor structure for docking purposes [36,37].ADFR suite [38] and Meeko python package were used to prepare the protein and ligand structures.AutoDock Vina v1.2.0 [39][40][41] was used in the docking runs.Details of docking setups can be found in Supplementary Table S4.Subsequently, after Tier I screening, we ranked compounds based on theoretical binding affinity (t-BA) scores and chose the top 2000 hit compounds for Tier II screening.
In Tier II screening, the top 2000 hits were further docked targeting SHP2 closed states with the SHP099 inhibitor, i.e., wtSHP2-b and mtSHP2-b.Again, the most dominant MD solution conformation extracted from the k-mean clustering analysis for each complex was utilized as the receptor structure for docking purposes.Subsequently, compounds were ranked based on t-BA-scores (Supplementary Table S2A,B).
For consensus docking results, we implemented two distinct selection criteria.The first criterion (t-BA) relies on t-BA-scores to prioritize compounds.Here we sorted compounds based on t-BA-scores for both receptors and then selected common compounds in the top 20 lists.This leads to six common top-hit compounds.We continued by selecting common compounds in the top 40 lists that are not already selected in the top 20 lists.This leads to 11 additional common top-hit compounds.We continued this way for the top 60, top 80, and top 100 lists, leading to a total of 73 consensus top hits (Supplementary Table S3).We further visualized these top hits and selected 18 potential compounds for further analysis (Table 1).In the h-VI criterion, we examined the spatial orientation of compounds within the allo-closite pocket and their interactions with residues from all three domains of SHP2.We not only visualized the top 2000 compounds, but also lower-ranked compounds.Finally, we found a total of 21 potential compounds (Table 2).

MD Simulations
Four crystal structures were adopted for MD simulations, including ligand-free wild type (PDB ID: 2SHP [10]) and mutant proteins (PDB ID: 6CRF [11]), and ligand-bound wild type (PDBID: 5EHR [8]) and mutant proteins (PDB ID: 6CRG [11]).The missing loops of all the crystal structures were modeled using Modeller in Chimera [42][43][44][45].The ff19SB force field was used to model proteins [46] and the TIP3P model [47] was used for water.The Amber-compatible GAFF2 force field parameters [48] for all ligands were generated with the Antechamber program of the AmberTools software suite version 2023 [37].All systems were solvated in a truncated octahedral box with a buffer of 8.0 Å and a NaCl concentration of 150 mM in the LEaP program of AmberTools 2023 [37].A 9.0 Å cutoff was applied for the short-range electrostatic and nonbonded interactions.The long-range electrostatics was handled with the PME approach with default settings [49].All bonds to the hydrogen atoms are constrained with the SHAKE algorithm [50,51] so that the leap-frog method [52] with a 2 fs time step could be used for time integration.All simulations were conducted with the PMEMD program for either the CPU or the GPU platform from the Amber software suite version 2023 [36,53].
For each simulation, the initial minimization involved 2000 steps via the steepest descent method, followed by an equal number of steps using the conjugate gradient technique.During minimization, all protein and ligand heavy atoms were restrained with a force constant of 10 kcal/mol-Å while hydrogen atoms and the solvents, including both water molecules and ions, are allowed to relax.Next, a heating step over 200 ps was conducted, starting from 0 K to 298 K with the Langevin thermostat [54] with a collision frequency of 1.0 per ps.Here the constant volume ensemble was used.Following the heating phase, an equilibration phase was undertaken for 5 ns under the NPT ensemble with the Berendsen barostat with default parameters [55], with the restraining force constant reduced to 5 kcal/mol-Å.The equilibration process was then repeated without any restraint.This was followed by the final production run up to 500 ns in the NVT ensemble in the Berendsen thermostat [55].MD simulations were repeated up to three times to collect sufficient amount of data for conformational analysis.All MD trajectory analysis was conducted with the CPPTRAJ program [56].Specifically, the k-means clustering method was used to extract dominant conformations from all trajectories and to detect conformational changes.Clustering number was set to 10.The heavy atom (up to side-chain CB atoms) RMSD was used as the distance metric.

Conclusions and Future Directions
In this study, molecular dynamics (MD) simulations were utilized to elucidate the structural and dynamic properties of SHP2 in both the wild-type (wtSHP2, closed state) and the E76K mutant (mtSHP2, open state) conformations, examining monomeric and dimeric forms as depicted in crystal structures.A comprehensive dataset was amassed from 28.6-µs MD simulations in an explicit solvent model across various systems.
Our findings reveal that each domain of wtSHP2 retains substantial stability, whether in monomer or dimer configuration.This stability is postulated to derive from the dense organization of wtSHP2, where domains are tightly packed against each other.Transitioning to mtSHP2, the N-SH2 and PTP domains maintain stability comparable to wtSHP2.However, the C-SH2 domain demonstrates increased flexibility, consistent with its less defined electron density in the crystal structure, suggesting looser packing both intramolecularly and intermolecularly.Notably, the crystal packing within the unit cell is significantly different between wtSHP2 and mtSHP2.In wtSHP2, our MD simulations of solvated monomers/dimers correspond well with the crystal structure.In contrast, the mtSHP2 presents a more complex packing, which may account for the disparities in domain stability observed in our simulations.
Building on these structural insights, our simulations highlight the impacts of the E76K mutation on the second allosteric binding site, where we observed a reduced binding efficacy due to the site's flat topology, which likely leads to increased inhibitor dissociation.A pivotal aspect of our research was the identification of a new secondary binding pocket, with the hope for robust and comprehensive inhibition of SHP2, especially in its mutant form.This allo-closite pocket, situated at the interfaces of the N-SH2, C-SH2, and PTP domains and remote from the SHP099 binding site, exhibited enhanced stability upon SHP099 binding in comparison to the apo form of wtSHP2, signifying the pocket's potential as a selective target that operates via a conformation-selection mechanism.
To further substantiate our findings, we conducted a virtual screening of the predicted pocket against the NCI cancer-specific small molecule library, implementing a three-stage screening protocol.The screening process yielded 39 promising compounds, identified through a consensus selection method complemented by human visualization after the initial screening phases.We then performed 100 ns MD simulations to weed out weak binders among the selected compounds within the allo-closite pocket.From these simulations, 11 compounds emerged as highly stable, signifying their potential as hit compounds.
Looking forward, the intricate structural understanding of SHP2 and its mutants obtained from this study opens new avenues for further research into novel therapeutic interventions.Our findings underscore the potential of allosteric inhibitors, particularly highlighting the importance of prolonged target binding for effective SHP2 suppression.While our work is purely computational and requires experimental validation, the identification and virtual screening validation of a novel allosteric pocket provides a compelling hypothesis for future investigation.We encourage experimentalists within the scientific community to explore this avenue further, following established strategies [11,12].Continued research in this direction promises to refine drug design strategies, potentially leading to more effective treatments for SHP2-associated cancers.

Figure 1 .
Figure 1.Overview of Structural Dynamics of the SHP2 Protein.(A) depicts the domain organization of the full-length SHP2 protein.(B) illustrates SHP2 in its closed or basal state, PDB ID: 2SHP [10].In the closed state, the DE-loop of the N-SH2 domain occludes the PTP domain catalytic cleft, effectively preventing the binding of diverse phosphorylated substrates and small molecules, and affecting the SHP2 function.(C) showcases the open state of SHP2, PDB ID: 6CRF [11], induced by the oncogenic mutation (E76K).The mutation causes a 120-degree C-SH2 domain rotation, exposing the catalytic cleft and shifting the N-SH2 domain across the PTP domain.In the open state, the DEloop of the N-SH2 domain no longer interacts with the catalytic cleft of the PTP domain, resulting in complete access to the catalytic cleft for pY-ligands.Additionally, this gain-of-function mutation compromises the binding cavity for SHP099.

Figure 1 .
Figure 1.Overview of Structural Dynamics of the SHP2 Protein.(A) depicts the domain organization of the full-length SHP2 protein.(B) illustrates SHP2 in its closed or basal state, PDB ID: 2SHP [10].In the closed state, the DE-loop of the N-SH2 domain occludes the PTP domain catalytic cleft, effectively preventing the binding of diverse phosphorylated substrates and small molecules, and affecting the SHP2 function.(C) showcases the open state of SHP2, PDB ID: 6CRF [11], induced by the oncogenic mutation (E76K).The mutation causes a 120-degree C-SH2 domain rotation, exposing the catalytic cleft and shifting the N-SH2 domain across the PTP domain.In the open state, the DE-loop of the N-SH2 domain no longer interacts with the catalytic cleft of the PTP domain, resulting in complete access to the catalytic cleft for pY-ligands.Additionally, this gain-of-function mutation compromises the binding cavity for SHP099.

Figure 2 .
Figure2.Analysis of SHP099 binding site dynamics across all SHP2 simulations.(A,B) depict SHP2 in its closed (wtSHP2) and open (mtSHP2) conformations.In the closed state, SHP099 binds from the back side, while the front-side access is hindered by a short disordered loop.However, mutation E76K alters the binding pocket through a rotation in the C-SH2 domain, resulting in the displacement of key residues Arg111 and Phe113, which play a role in the creation of the binding pocket for SHP2 inhibitors.(C) represents the CA distances between Arg111 and Phe113 that are involved in SHP099 binding across all SHP2 configurations.(D,E) illustrate the superimposed structures from the top MD solution clusters of each SHP2 system, while the binding site residues form a consistent cavity in the wtSHP2 and mtSHP2-b states, a discernible shift in two residues from the C-SH2 domain away from the binding pocket is evident for the pocket disruption in the mtSHP2-ub conformation.The green sticks denote the crystallographic structure from the wtSHP2-b state.

Figure 2 .
Figure 2. Analysis of SHP099 binding site dynamics across all SHP2 simulations.(A,B) depict SHP2 in its closed (wtSHP2) and open (mtSHP2) conformations.In the closed state, SHP099 binds from the back side, while the front-side access is hindered by a short disordered loop.However, mutation E76K alters the binding pocket through a rotation in the C-SH2 domain, resulting in the displacement of key residues Arg111 and Phe113, which play a role in the creation of the binding pocket for SHP2 inhibitors.(C) represents the CA distances between Arg111 and Phe113 that are involved in SHP099 binding across all SHP2 configurations.(D,E)illustrate the superimposed structures from the top MD solution clusters of each SHP2 system, while the binding site residues form a consistent cavity in the wtSHP2 and mtSHP2-b states, a discernible shift in two residues from the C-SH2 domain away from the binding pocket is evident for the pocket disruption in the mtSHP2-ub conformation.The green sticks denote the crystallographic structure from the wtSHP2-b state.

Figure 3 .
Figure 3. Conformational Dynamics of the Allosteric Site-2 (allo-site-2) Pocket in SHP2.(A) shows the superimposition of all closed-state MD solution structures, including wtSHP2-ub/b and mtSHP2b, along with their corresponding crystal structures 2SHP, 5EHR, and 6CRG, respectively, revealing the pocket's conformational consistency.(B) provides a distance analysis between two salt-bridging residues flanking pocket allo-site-2.In wtSHP2-b, the two residues start from their initial stable state but transition into a disrupted state, leading to a flat binding surface.In mtSHP2-b, the two residues transition into a closed state, forming a salt bridge and capping the binding pocket.(C,D) illustrate conformational variations observed in the simulations: Shift 1 in a wtSHP2-b simulation and Shift 2 in an mtSHP2-b simulation.(E) showcases the surface representation of allo-site-2.In MD solution structures, the disordered loop stays either in an 'in' conformation in wt-SHP2-b, capping allo-site-2, or in an 'out' conformation in mt-SHP2-b, allowing ligand-binding.The three crystal structures used in MD simulations are shown along with the two MD solution structures.Note that the disordered loop is missing in all three crystal structures.

Figure 3 .
Figure 3. Conformational Dynamics of the Allosteric Site-2 (allo-site-2) Pocket in SHP2.(A) shows the superimposition of all closed-state MD solution structures, including wtSHP2-ub/b and mtSHP2-b, along with their corresponding crystal structures 2SHP, 5EHR, and 6CRG, respectively, revealing the pocket's conformational consistency.(B) provides a distance analysis between two salt-bridging residues flanking pocket allo-site-2.In wtSHP2-b, the two residues start from their initial stable state but transition into a disrupted state, leading to a flat binding surface.In mtSHP2-b, the two residues transition into a closed state, forming a salt bridge and capping the binding pocket.(C,D) illustrate conformational variations observed in the simulations: Shift 1 in a wtSHP2-b simulation and Shift 2 in an mtSHP2-b simulation.(E) showcases the surface representation of allo-site-2.In MD solution structures, the disordered loop stays either in an 'in' conformation in wt-SHP2-b, capping allo-site-2, or in an 'out' conformation in mt-SHP2-b, allowing ligand-binding.The three crystal structures used in MD simulations are shown along with the two MD solution structures.Note that the disordered loop is missing in all three crystal structures.

Figure 4 .
Figure 4. Dynamics of the allo-closite Pocket in SHP2 States.(A) illustrates the Cα distance (CA-dist) among key residues from the SH2 and PTP domains that constitute the allo-closite pocket in SHP2 simulations unbound to SHP099.Variations in CA-dist provide insights into the pocket's integrity, where increased distances suggest disruption, and decreased distances indicate pocket formation.Converging CA-dist values are indicative of allo-closite pocket establishment.(B) displays superimposed structural snapshots demonstrating the allo-closite pocket's conformational shifts in SHP099free SHP2 simulations.(C) showcases superimposed snapshots in SHP99-bound SHP2 simulations, underscoring the enhanced stability and potential therapeutic significance of the allo-closite pocket.(D) presents the CA distance (CA-dist) plot for SHP2 simulations when bound to SHP099, which aids in assessing the pocket dynamics upon inhibitor binding.

Figure 4 .
Figure 4. Dynamics of the allo-closite Pocket in SHP2 States.(A) illustrates the Cα distance (CAdist) among key residues from the SH2 and PTP domains that constitute the allo-closite pocket in SHP2 simulations unbound to SHP099.Variations in CA-dist provide insights into the pocket's integrity, where increased distances suggest disruption, and decreased distances indicate pocket formation.Converging CA-dist values are indicative of allo-closite pocket establishment.(B) displays superimposed structural snapshots demonstrating the allo-closite pocket's conformational shifts in SHP099-free SHP2 simulations.(C) showcases superimposed snapshots in SHP99-bound SHP2 simulations, underscoring the enhanced stability and potential therapeutic significance of the alloclosite pocket.(D) presents the CA distance (CA-dist) plot for SHP2 simulations when bound to SHP099, which aids in assessing the pocket dynamics upon inhibitor binding.

Figure 5 .
Figure 5. Tier III MD screening of common compounds for SHP2 closed states.(A) Ligand R for wtSHP2 and mtSHP2 based on t-BA criteria.(B) Ligand RMSD based on h-VI criteria.(C) Sta distributions in collective MDs.(D) Superposition of all ligands to show the ligand distributio the allo-closite pocket.(E) Conformational change of the allo-closite gate-keeper residues, for salt bridge upon binding of hit compounds.(F) 2D structures of 11 stable common compound

Figure 5 .
Figure 5. Tier III MD screening of common compounds for SHP2 closed states.(A) Ligand RMSD for wtSHP2 and mtSHP2 based on t-BA criteria.(B) Ligand RMSD based on h-VI criteria.(C) Stability distributions in collective MDs.(D) Superposition of all ligands to show the ligand distributions in the allo-closite pocket.(E) Conformational change of the allo-closite gate-keeper residues, forming salt bridge upon binding of hit compounds.(F) 2D structures of 11 stable common compounds.

Author Contributions:
Conceptualization, A.U.R. and R.L.; data collection, Y.W. and C.Z.; validation, C.Z., Q.Z. and R.L.; visualization and data analysis, A.U.R. and C.Z.; writing-original draft preparation, A.U.R. and C.Z.; writing-review and editing, A.U.R., C.Z., Y.W., Q.Z. and R.L; project administration, R.L.; funding acquisition, R.L.All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by NIH/NIGMS grant number GM130367 to R.L. Data Availability Statement: Simulation data available from authors upon request.

Table 1 .
The top 18 compounds selected from 73 common top-hit compounds at the end of Tier II screening (affinity in kcal/mol).

Table 1 .
The top 18 compounds selected from 73 common top-hit compounds at the end of Tier II screening (affinity in kcal/mol).

Table 2 .
Manually selected top 21 compounds from the Tier I and Tier II screening (affinity in kcal/mol).