SHP099

Probing the acting mode and advantages of RMC-4550 as a src-homology 2 domain-containing protein tyrosine phosphatase (SHP2) inhibitor at molecular level through molecular docking and molecular dynamics

Abstract
The over-activation of Ras/mitogen activated protein kinase (MAPK) signaling pathway associated with a variety of cancers is usually related with abnormal activation of Src-homology 2 domain-containing protein tyrosine phosphatase (SHP2). For this purpose, SHP2 has attracted extensive interest as a potential target for cancer treatment. RMC-4550, as a newly developed selective inhibitor of SHP2, possesses of an overwhelming advantage over the previous generation inhibitor SHP099 in terms of vitro activity. However, the binding mode of SHP2 with RMC-4550 and the reason for the high efficiency of RMC-4550 as SHP2 inhibitor at molecular level are still unclear. Therefore, in this study, the binding mode of RMC-4550 with SHP2 and the superiorities of RMC-4550 as inhibitor at binding affinity and dynamic interactive behavior with SHP2 were probed by molecular docking and molecular dynamics simulations. By comparing the results of molecular docking, it was found that SHP2 formed more tight interaction with RMC-4550 than that with SHP099. Subsequently, a series of post-dynamic analyses on three simulation trajectories (SHP2WT, SHP2SHP099 and SHP2RMC-4550) were performed and found that the SHP2 protein bound with RMC-4550 maintained a firmer interaction between N-Src-homology 2 (N-SH2) and PTP-domain throughout the MD simulation, leading to a more stable protein conformation. The finding here provides new clues for the design of SHP2 inhibitor against the over-activation of Ras/MAPK pathway.CHARMM (Chemistry at HARvard Macro-molecular Mechanics); DCCM (dynamic cross-correlation mapping); DSSP (definition of secondary structure of proteins); GEFs (guanine nucleotide exchange factors); GAPs (GTPase-activating proteins); MAPK (Mitogen activated protein kinase); MD (Molecular dynamics); MM-PBSA (Molecular Mechanics Poisson Boltzmann Surface Area); PTPN11 (protein tyrosine phosphatase 11); PCA (principal component analysis); PDB (protein data bank); RIN (residue interaction network); RING (Residue Interaction Network Generator); RMSD (Root mean square deviation); RMSF (Root mean square fluctuation); RTKs (Receptor tyrosine kinases); SH2 (Src homology 2); SHP2 (Src-homology 2 domain-containing protein tyrosine phosphatase); SPC(single-point charge)

Introduction
Mitogen activated protein kinase (MAPK) pathway is a vital signal transduction system for eukaryotic cell to regulate intracellular responses on the basis of extracellular signal(Easton, Royer, & Middlemas, 2006). The activation of MAPK relies on the interrupt of the Ras guanine nucleotide exchange cycle regulated by guanine nucleotide exchange factors (GEFs) and GTPase-activating proteins (GAPs)(Heuberger et al., 2014). However, the abnormal activation of Ras/MAPK signaling pathway can trigger the excessive division of cell. Therefore, the aberrant signaling of Ras/MAPK has been defined as the main mechanism in the many malignant tumor, including mammary cancer, gastric cancer, melanoma, oophoroma, juvenile myelomonocytic leukemia(Y. Patel et al., 2016). One of the common causes of tumor-driven hyperactivation of Ras/MAPK pathway is abnormal activation in downstream component of the pathway such as receptor tyrosine kinases (RTKs)(Hu et al., 2018).Remarkably, the hypernomic level of phosphorylation in Src-homology 2 (SH2) domain-containing protein tyrosine phosphatase (SHP2) protein has been detected in many patients with diseases related to Ras/MAPK pathway(S. J. Yu et al., 2011). As a signal-transducer downstream of multiple RTKs, SHP2 can promote Ras/MAPK activation by integrating growth factor signals. SHP2, encoded by protein tyrosine phosphatase 11 (PTPN11), comprises two Src-homology 2 domains (N-SH2 domain and C-SH2 domain), a catalytic PTP domain and a C-terminal tail with the binding site for tyrosyl phosphorylation (C. Chen et al., 2015; Fragale, Tartaglia, Wu, & Gelb, 2004). Due to the intramolecular interaction between the N-SH2 domain and PTP domain, SHP2 is in a self-inhibition state with low catalytic activity. Once the phosphotyrosyl-containing proteins bind to the SHP2 at separate site on the N-SH2 and C-SH2 domains, the autoinhibitory interaction of the SHP2 is disrupted(J. Xie et al., 2017). Then the changed conformation can release the N-SH2 domain from the PTP domain and the sites with catalytic active located in the PTP domain will be exposed to the substrate more easily, which can lead to the increased activity of the protein (Scheme 1)(Butterworth, Overduin, & Barr, 2014; Oishi et al., 2009). Therefore, keeping the low activity of protein by maintaining the closed structure of protein is the key factor for the design of SHP2 inhibitors.

Along with the computer-aided technology widely being used in receptor-based drug designing field, the molecular docking is considered as a reliable tool in predicting the binding mode and affinity between protein and ligand(Malathi, Ramaiah, & Anbarasu, 2019). By binding the ligand to a known active site, the most suitable conformation of ligand is determined according to orientation and interaction(J. Wang et al., 2016). Moreover, in recent years, molecular dynamics (MD) simulation has advanced into a mature technology for exploring the complicated biologic problems that are difficult to be disposed of experimentally(Hamed & Arya, 2016; Mahmoodi, Mehrnejad, & Khalifeh, 2018). Based on numeric integration of Newtonian motion formulas, MD simulation attempts to gain the energetically permissible conformational state of a protein or a protein–ligand complex and provides us extremely credible information to evaluate the effects of individual ligand on protein receptor during entire simulation(Hatmal, Jaber, & Taha, 2016). Therefore, in this study, the molecular docking serves as a tool to probe the differences in mode and affinity of interacting with SHP2 between RMC-4550 and SHP099. Meanwhile, the effects of ligand upon protein conformation and residue-residue interaction can be monitored by 150 ns MD simulation.After the molecular docking being completed, by constructing 2D diagrams of protein-ligand interaction, the docking results were visualized to compare the binding modes of the SHP2 with SHP099 and RMC-4550 ligands.

And the binding free energy between individual residue and ligand was calculated by Molecular Mechanics Poisson Boltzmann solvent-accessible Surface Area (MM/PBSA) to predict the affinity of interaction between SHP2 and two ligands, respectively. Subsequently, a series of post-dynamic analyses were performed on 150 ns trajectories of wild type (SHP2WT) system, complex of SHP2 with ligand SHP099 (SHP2SHP099) system and complex of SHP2 with ligand RMC-4550 (SHP2RMC-4550) system, respectively. Through the studies on stability, flexibility, conformational state, correlation motion of residue-residue, secondary structure and residue interaction network, the effects of two ligands on the dynamics of the protein were also understood. Finally, the information from molecular docking and dynamic simulations were gathered to compare the difference in mode and affinity of interaction with SHP2 between RMC-4550 and SHP099, and to clarify the superiority of RMC-4550 over SHP099 at the molecular level. The results from our study can give valuable insights into the allosteric inhibitory mechanism of SHP2 protein and provide the useful clues for the design of SHP2 inhibitor against Ras/MAPK pathway.The high‑ resolution crystallography structure of SHP2 protein was retrieved from protein data bank (identifier PDB ID: 5EHR), with resolution and residue count of 1.7 Å and 526, respectively. To optimize the initial target, the crystal structure was processed by the clean protein protocol of Protein Preparation workflow to repair and refined the missed residues and atoms, delete the nonprotein atoms (including water molecules and ligand), assign bond orders, add the hydrogen and remove chain B. Subsequently, the ligands were fully prepared by dislodging duplicates, enumerating tautomer and stereo isomers, maintaining ionization state and generating appropriate 3D conformation. In this study, all the preparations above were completed by Discovery Studio v3.5.

Molecular docking processes were completed by the LigandFit protocol in Discovery Studio v3.5.First, a sphere with a volume of 10 Å was generated based on the binding site of co-crystallized ligand in the used crystal structure. And then, the different conformations were forwarded onto the binding site. Therewith, the docking poses were minimized with Chemistry at HARvard Macro-molecular Mechanics (CHARMM) forcefield(S. Patel, Mackerell, & Brooks, 2004). After docking, the conformations were filtered according to the Glide Score value and the matched degree in orientation between the docked ligand conformation and conformation in crystal complex(K. C. Chen, Sun, & Chen, 2014; Sun, Zheng, & Zhang, 2017). The most optimal conformation was selected for analysis of binding pose, and the detailed interactions (hydrogen bonds and hydrophobic) between protein and ligands were shown as 2D diagrams constructed by Protein and Ligand module in Discovery Studio v3.5.MD simulations for SHP2WT, SHP2SHP099 and SHP2RMC-4550 systems were performed by using the ‘GROMACS 4.5.5’ software package with GROMACS 43a1 force field(Pol-Fachin, Fernandes, & Verli, 2009). Each system was simulated in a closed dodecahedral box and dissolved with single-point charge (SPC) water molecules.(Takano, Nakata, Yonezawa, & Nakamura, 2016). Herein, the minimum distance from the solute to the edge of the box in the system was limited to1.0 nm.

For neutralizing the charges, the appropriate number of counter ions was added to the system. Subsequently, the steepest descent method was used for minimizing the energy of the system to 1000 kJ·mol-1/nm to ensure that there was no steric collision and inappropriate geometry in the initial structure(L. Wang et al., 2018). Afterwards, under the 100 ps NVT and NPT simulation condition, the temperature and pressure of systems were stable at 300 K and 1 br, respectively(Bhutani, Loharch, Gupta, Madathil, & Parkesh, 2015). Finally, each system was simulated for 150 ns, and each trajectory was sampled every 20 ps interval.Principal component analysis (PCA), as a statistical method, is utilized for obtaining concerted movement of proteins(Karunakaran & Srikumar, 2018). This method is widely used to reduce the dimensionality of the data gained from molecular dynamics simulations and to screen observed motions from the largest to smallest spatial scales through the decomposition process(David & Jacobs, 2014; Kan, Fang, Chen, Wang, & Deng, 2016). Based on the ensemble formula, PCA can be used to gain a covariance matrix with elements Cij of coordinates i and j, and the formula is given as:Cij =〈( xi -〈xi〉) (xj -〈xj〉)〉(i, j = 1, 2, 3,…,N) Here, xi represents the Cartesian coordinates of the ith Ca atom and xj represents the Cartesian coordinates of the jth Ca atom; and the Angular brackets “〈〉” denote the time average of all relevant configurations during the entire simulation; and i and j means the number of Ca atom(Yesudhas et al., 2016; Zhou, Zhang, Chen, Zhao, & Zhong, 2016). The dynamic motion of the atoms in the system is calculated from the simulation trajectories to analyze the conformational differences between the system during the simulation(Fakhar et al., 2017; Yesudhas et al., 2016).

Results and discussion
The docking results of SHP2WT with RMC-4550 and SHP099 were completed in the Discovery Studio v3.5, respectively, and the corresponding docking poses were visualized in workspace. The top-ranked docking score for RMC-4550 and SHP099 was 46.6943 and 40.9327 respectively, indicating that RMC-4550 possessed of more appropriate molecular orientation for the binding site in SHP2(Luo et al., 2012). Subsequently, the binding pockets generated by the interaction surface area were displayed in Figure 1. From the Figure 1, it was revealed that both ligands were embedded between the C-SH2 (residues Leu103-Lru216) and the PTP domain (residues Asn217- Thr524), therefore they were considered for sharing the same binding site. The result above suggested a similar action mechanism for RMC-4550 and SHP099 as SHP2 inhibitors, and namely SHP2 protein might be maintained in a stability “closed” structure by locking the C-SH2 and PTP domain. electrostatic bond interactions. Moreover, more residues in SHP2 were involved in forming the interaction with RMC-4550 than that with SHP099. Within the docking binding site of SHP2, compared to SHP099, additional hydrogen bonds were formed between the hydroxyl provided by the RMC-4550 and the residues Thr218 and Thr219, which might explain the advantage of introducing a hydroxyl on ortho-position to replace the amino. In addition, residues Tpr112, Pro215 and Tyr511 could provide the extra Van der Waals interactions to stabilize the SHP2RMC-4550 complex. Meanwhile, the existence of a charged polar residue Lru216 could play a positive role in stabilizing the SHP2RMC-4550 complex via electrostatic interaction. The above analysis demonstrated that more residues in SHP2 provided interactions with RMC-4550, which might lead to a more stable binding of RMC-4550 to SHP2 and thus RMC-4550 could better maintain the “closed” structure of SHP2 to the SHP2, and herein, Van der Waals energies were the main forces for driving the ligands binding to SHP2. Subsequently, to figure out the different component contributions to the interactions of SHP2 with ligands, the binding free energies were decomposed into the per-residue coming from the residues side-chain and backbone to obtain the differences in source of binding affinity between SHP2 and two ligands(Tang & Chen, 2015; Zhang et al., 2017). The results from Figure 3 demonstrated that, in SHP2SHP099 complex, the most key contributions to the binding came from the residues His114 (from C-SH2), residues Thr218, Glu249, Gln257, Pro491, Lys492 and Gln495 (from PTP domain). Whereas, in SHP2RMC-4550 complex, the residues that were favorable for binding were Glu110, Trp112, Pro215, Asn217, Glu250, Gln257, Asp489, Pro491, Lys492 and Tyr511. Consequently, the increased van der Waals interactions were beneficial to stabilize the binding of SHP2 and RMC-4550, which might be one of the reasons why RMC-4550 had better effect of inhibition than SHP099.

Hydrogen bond was a ubiquitous chemical interaction whose formation played an enormouse role in stabilizing protein-ligand complexes(Sigala et al., 2015). The comprehensive investigation into occupancy and distance of hydrogen bond was thus necessary to further understand the contribution of per-residue to the affinity of ligand with receptor. In this work, the distance between the donor and acceptor atoms forming the hydrogen bond was recorded, and the hydrogen bond occupancy was calculated based on the ratio of hydrogen bond formed during the entire MD simulation. The cutoff of hydrogen bond distance was set between 0.25 nm and 0.35 nm and all the atoms that could form hydrogen bond were taken into consideration(H. Xie et al., 2015). And then, all the residues participating in the interactions of hydrogen bond were displayed in Figure 4. From the Figure 4A, the residues that formed hydrogen bond with SHP099 were Glu110, Arg111, Phe113, Glu250. While, from the Figure 4B, the RMC-4550 mainly formed interaction of hydrogen bond with residues Arg111, Phe113, His114, Thr218, Thr219. It was obvious that the two residues Arg111 and Phe113 forming hydrogen bonds in both complexe played a significant role in interaction between ligands and allosteric site of SHP2. It’s worth noting that the additional hydrogen bonds formed with residues His114, Thr218, Thr219 might be another reason for higher inhibitory effect of RMC-4550. The hydrogen bond analyzed above maintained over 80% occupancy throughout the MD simulation, elucidating that the hydrogen bonds were stably formed in complexes, which was in accordance with the Figure 2. The result from analyses on hydrogen bond revealed that the two inhibitors existed the different binding mode and the hydrogen bond have obviouse effect on the binding affinity.

The crystal structure of SHP2 protein from x-ray demonstrated that, due to the interaction between N-SH2 and PTP domain, the protein was in an autoinhibited state with low activity. SHP2 bound with Ligands at the allosteric site could stabilize the autoinhibited (closed) conformation of the protein (Xue, Yang, Wang, Liu, & Yao, 2014). To evaluate the effects of two ligands on protein stability, the root mean square distance (RMSD) of the backbone atoms relative to initial structure over the course of the simulations was monitored and record in Figure 5(Timson & Lindert, 2013). As illustrated in Figure 5, the values of RMSD for SHP2WT, SHP2SHP099 and SHP2RMC-4550 systems tended to be convergent after ∼10 ns, indicating that the three simulated systems attained sufficient equilibrium in the last 140 ns of MD simulations. Therefore, the trajectories of 10-150 ns were intercepted for further post dynamic analyses. Subsequently, by calculating the average RMSD value of each system, it was observed that the SHP2WT system with a greater value of RMSD curve underwent more largely structural alteration, while the systems bound with two inhibitors remained the relatively stable conformation during the simulation. Moreover, the lowest RMSD value in SHP2RMC-4550 complex suggested that RMC-4550 was provided with a greater advantage in stabilizing the conformation of the protein than SHP099 system. This revealed that the binding of SHP2 with ligands promoted the stabilization of the crucial interactions between N-terminal and PTP domain, and consequently increasing the rigidity of these loops. Moreover, compared to the SHP2SHP099 system, the SHP2RMC-4550 system had lower fluctuations of domains, which demonstrated that RMC-4550 could better match allosteric site of SHP2 and achieved stronger inhibitory effect projections of first two major eigenvalues (PC1 and PC2) for the three systems were evident in Figure 7. The result from Figure 7 highlighted that all three systems converged from energetically unstable conformation state (shown in scattered blue dots) to a stable conformation state (displayed in compact red dots).

Consequently, in the process of simulations, these systems experienced non-periodic jump between these two conformational states. The further observation on the PCA scatter plots revealed the conformational distribution of the SHP2WT system along the first two principal components was irregular and disorder (Figure 7A), whereas the conformational distributions of the SHP2SHP099 and SHP2RMC-4550 systems were evenly distributed on both sides of the diagonal (Figure 7B and Figure 7C), which indicated that the motion of SHP2 was restricted after bound with ligands. Thereafter, the distinct differences in eigenvalue calculated from 10-150 ns MD trajectories were observed between SHP2WT, SHP2SHP099 and SHP2RMC-4550 systems, indicating that there was significant difference in the motions of these three systems. Compared to SHP2WT and SHP2SHP099 systems, SHP2RMC-4550 system was found to display the highest correlated motion along PC1 and PC2 with the variance proportion of 46.5% and 14.7%, respectively. This suggested that SHP2 bound with RMC-4550 had a lesser occupancy of phase space. The results above might be caused by more compact structure and larger molecule rigidity of the SHP2 after bound with RMC-4550, which provided a solid evidence for the superiority of RMC-4550 as an inhibitor over SHP099 the correlated motion between residue pairs(C. Chen et al., 2018). Herein, the positive correlation (from 0 to 1) meant that the residues moved in the same direction while the negative correlation (from -1 to 0) represented that the residues moved in the opposite direction(R. R. Wang et al., 2019).

The positive correlation and the negative correlation were exhibited as green areas and red areas on Figure 8, respectively. As shown in Figure 8, after SHP2 bound with ligands, the D’E-loop showed the enhanced positive correlation with pTyr-Loop, P-loop and Q-loop, and the weakened negative correlation with WPD-loop. This suggested that the interaction between the N-SH2 and the PTP-domain was more stable in the SHP2SHP099 and SHP2RMC-4550 complexes, allowing the protein to maintain a closed structure. Furthermore, it was found from further observation that the D’E-loop displayed stronger positive-correlated motion with PTP-domain in SHP2RMC-4550 complex than that in SHP2SHP099 complex, which indicated that RMC-4550 had a greater advantage in stabilizing the conformation of SHP2 protein by fastening the interaction between N-SH2 and PTP-domain than SHP099 secondary structure during entire simulation was examined by Dictionary of Secondary Structure of Proteins (DSSP) tool(Rashid, Saraswathi, Kloczkowski, Sundaram, & Kolinski, 2016). The elements of secondary structure in each system were shown as Figure 9. The Figure S1 obtained from PDBsum showed that the binding of the ligands did not alter the secondary structure type of the whole protein. However, it was worth noting from Figure 9 that the secondary structure of the key regions tended to remain stable after bound with the ligands.The results from Figure 10 showed that the structures of SHP2SHP099 and SHP2RMC-4550 complexes were more compact than that of SHP2WT system. Subsequently, in order to facilitate observation, the region where the N-SH2 directly interacted with the PTP-domain was amplified and circled with the green elliptic frame. As could be seen from the figure 10, residues located at N-SH2 (marked in red) interacted with PTP region to achieve low activity of self-inhibition, which was consistent with previous report. And then, by comparing RINs of three systems, the interactions between the key regions in the SHP2SHP099 and SHP2RMC-4550 systems were increased significantly due to the existences of ligands. This suggests that the ligands could maintain the autoinhibited conformation of SHP2 by firming the interaction between N-SH2 and PTP domain, which was the detailed reason for inhibitory effect of SHP2 allosteric inhibitors. Moreover, more interactions in residues Asn58, Thr59, Asp61, Tyr63, Ala72, Leu74 were observed in SHP2RMC-4550 complex compared to SHP2SHP099. This illustrated that RMC-4550 held an obvious advantage over SHP099 in stabilizing internal self-interaction of SHP2 protein. Therefore, the result from RINs might clarify the most detailed reason at molecular level for the high activity of RMC-4550 as SHP2 inhibitor.

Conclusion
This paper mainly emphasized on probing the binding mode of known ligand RMC-4550 with SHP2 and the advantages of RMC-4550 as inhibitor at molecular level by molecular docking and molecular dynamics simulation. The result from molecular docking showed that RMC-4550 was bound well to the binding pocket of SHP2 by mainly forming interactions of hydrogen bond and Van der Waal. Besides, RMC-4550 interacted with more residues of SHP2, indicating that SHP2 bound with RMC-4550 was more stable than that with SHP099. With the aid of calculating binding free energy and hydrogen bond distance, the stability of RMC-4550 bound to SHP2 was further verified. Thereafter, the analyses on the RMSD and RMSF indicated that, compared to SHP099, RMC-4550 could lead to the higher rigidity of protein structure and lower flexibility in key loop regions. This indicated that the SHP2 bound with RMC-4550 remained the original closed structure to the greatest extent during the entire dynamics simulation. And then, PCA and DCCM were conducted to compare the conformational states and correlative motions between residue pairs of the three systems, which found that the SHP2 bound with RMC-4550 was provided with the lowest occupancy of phase space, and the stronger positive correlations between N-SH2 and PTP domain were enhanced during simulation. The analyses above suggested that the interaction between N-SH2 and PTP domain within the SHP2 was more stable, leading to a more compact protein conformation. Finally, by constructing the RINs, it was found that, after SHP2 bound with RMC-4550, the D’E-loop penetrated into the PTP-domain was more stably due to the interaction with more residues, which might be the most detailed reason for the superiority of RMC-4550 as SHP2 inhibitor over SHP099. The finding of this paper explained the action mechanism and advantages of RMC-4550 as an allosteric inhibitor of SHP2 at the molecular level and provided a reliable clue for the design of new SHP099 inhibitor against Ras/MAPK pathway.