Characterizing the Hot Spots Involved in RON-MSPβ Complex Formation Using In Silico Alanine Scanning Mutagenesis and Molecular Dynamics Simulation

Purpose: Implication of protein-protein interactions (PPIs) in development of many diseases such as cancer makes them attractive for therapeutic intervention and rational drug design. RON (Recepteur d’Origine Nantais) tyrosine kinase receptor has gained considerable attention as promising target in cancer therapy. The activation of RON via its ligand, macrophage stimulation protein (MSP) is the most common mechanism of activation for this receptor. The aim of the current study was to perform in silico alanine scanning mutagenesis and to calculate binding energy for prediction of hot spots in protein-protein interface between RON and MSPβ chain (MSPβ). Methods: In this work the residues at the interface of RON-MSPβ complex were mutated to alanine and then molecular dynamics simulation was used to calculate binding free energy. Results: The results revealed that Gln193, Arg220, Glu287, Pro288, Glu289, and His424 residues from RON and Arg521, His528, Ser565, Glu658, and Arg683 from MSPβ may play important roles in protein-protein interaction between RON and MSP. Conclusion: Identification of these RON hot spots is important in designing anti-RON drugs when the aim is to disrupt RON-MSP interaction. In the same way, the acquired information regarding the critical amino acids of MSPβ can be used in the process of rational drug design for developing MSP antagonizing agents, the development of novel MSP mimicking peptides where inhibition of RON activation is required, and the design of experimental site directed mutagenesis studies.


Introduction
Protein-protein interactions (PPIs) are involved in many biological processes as key regulatory steps 1 and when aberrantly regulated are implicated in the development of many diseases such as cancer. 2-4 These versatile roles make PPIs attractive for therapeutic intervention and rational drug design. 5-7 Different classes of the approved therapeutic agents or those in development stages have been shown to interfere with PPIs in order to overcome the corresponding diseases. [8][9][10][11] In PPI, the amino acids at the interaction interface have great importance in terms of starting point for the initiation of the biological and cellular functions. 12 Usually several residues are exposed at the interface between the interacting proteins, but they do not contribute equally to the binding energy. The important key residues when mutated into alanine residue weaken the binding strength (increase of free energy of binding at least 2.0 kcal/mol) are called "hot spots". 13,14 Identification of the hot spots is critical in designing therapeutic agents which exert their effects by influencing PPIs. One of the experimental methods commonly used for identification of these hot spots is site-directed mutagenesis followed by comparative functional assays of the mutated and the wild type proteins; however, these experiments are timeconsuming and expensive. 15 Structural elucidation of the partner proteins within a complex by means of biophysical methods such as X-ray crystallography and NMR is also possible but again costly and demanding. 16 In the era of modern drug discovery and development, the use of in silico methods shortens the rational drug design process in terms of both time and cost. [17][18][19][20][21] In this regard, identifying hot spots is not an exception. 22,23 Computational alanine scanning mutagenesis is a virtual method which has been extensively used for the characterization and prediction of hot spots in proteinprotein, protein-DNA and protein-small molecule complexes. 22, [24][25][26][27][28][29] Charged, polar, or bulky amino acids are virtually mutated to a neutral, small and non-polar amino acid such as alanine and then binding free energy is calculated for both wild type and mutant forms in order to estimate the contribution of the mutated residues to the binding energy. 30,31 One of the most routinely used approaches for computational estimation of binding free energy is based on accessible surface area models of implicit solvation method where molecular mechanics data are treated by Generalized Born surface area (MM-GBSA) algorithm. [32][33][34][35][36][37][38][39][40] Tyrosine kinase receptors (TKRs) involved in well characterized protein-protein interactions are among potential candidate targets for anticancer drug development. [41][42][43][44][45][46] TKRs are cell surface receptors for different polypeptide ligands and have pivotal roles in regulation of many cellular functions and physiological events 47,48 and when aberrantly expressed and activated play key functions in development and progression of different types of cancers. [49][50][51][52][53][54] Ligand-mediated receptor dimerization is the main mechanism of activation triggered by ligand binding to the extracellular domain of its specific receptor. [55][56][57] This protein-protein interaction causes receptor dimerization followed by authophosphorylation of tyrosine residues located within the intracellular tyrosine kinase domain (catalytic tyrosines) followed by phosphorylation of tyrosine residues located within the C tail (docking tyrosines) that become the docking site for adaptor/effector proteins responsible of transducing the downstream signaling pathways resulting in cellular proliferation, differentiation, metabolism, survival, migration, and cell cycle control. 58 In principle, all PPIs mediated by TKRs (including the downstream PPIs) could be targeted for cancer therapy 59,60 but generally therapeutic PPI inhibitors interfere with the binding of endogenous ligands to the receptor. [61][62][63][64][65][66][67] Therefore, it is obvious that uncovering the details of PPIs between TKRs and their ligands can provide useful information applicable to design of new anticancer agents. RON (Recepteur d'Origine Nantais) is a member of TKRs superfamily, its role in tumorigenesis has been established in different cancer types and numerous studies have suggested RON as a promising target for anticancer drug development. 68,69 RON also known as MSTR1 (Macrophage Stimulating Receptor1) belongs to MET proto-oncogene family, 70 and is usually expressed at low levels in normal tissues while it is highly expressed in cancer cells. 71 Structurally, RON is a disulfide linked heterodimer protein made of two chains, an extracellular α-chain and a β-chain which consists extracellular, transmembrane, and intracellular regions. The extracellular domain comprises three distinct domains including Sema, Plexin-Semaphorin-Integrin (PSI), and three Immunoglobulin-Plexin-Transcription factor (IPT1-IPT3) domains. 68 The natural ligand of RON is MSP (Macrophage Stimulating Protein), 72 a member of plasminogen-related kringle protein family 73 which is a heterodimeric protein made of an α-chain composed of four kringle domains and a β-chain containing a serine protease-like domain. 74 The α-and βchains of MSP show low and high affinities to RON Sema domain, respectively. 75 Several monoclonal antibodies against RON extracellular domain have been developed (in preclinical phases) to specifically inhibit the protein-protein interactions between RON and MSP. 69 Identifying the key residues working as hot spots responsible for receptor-ligand (RON-MSP) interaction is of great importance for drug design and development.
The aim of the current study is to identify hot spots involved in RON-MSPβ interaction using in silico alanine scanning mutagenesis by MM-GBSA method. The results can be used in anticancer drug designing where inhibition of RON is needed.

Structure preparation and in silico alanine mutagenesis
Experimental coordinates of RON complexed with MSP (PDB ID: 4QT8) determined at 3.0 Å resolution by X-ray crystallography 76 was retrieved from the Protein Data Bank at the Research Collaboratory for Structural Bioinformatics (http://www.rcsb.org/pdb/ home/home.do). 77 Preparation of structures along with mutation of the residues were carried out using Swiss-Pdb Viewer (DeepView) version 4.01. 78 Only one of the complexes in the reported crystal structure was used (chains B and D) for further analysis. The residues at the RON-MSPβ interface were inferred based on crystal structure reported by Chao and collaborators 76 in both ligand and receptor were virtually mutated to alanine as listed in Table 1.

Ligand-receptor binding free energy calculations using MM-GBSA method
Energy minimization and binding free energy calculation were performed using the Assisted Model Building with Energy Refinement (AMBER) suite of programs (version 14) 79 operating on a Linux-based (Centus 6.8) GPU work station consisting of four Nvidia K40 M (each has 12 GB RAM and 2880 cuda cores), 2X Intel Xeon E5-2697 v2, 2.7 GHz (total of 48 cores), total RAM = 128 GB. The energy minimization of RON-MSP complex was carried out using AMBER-ff99 force field. 80 Briefly, the usable coordinate files for AMBER (i.e. *.prmtop and *.inpcrd) were generated using leap module. Then, a correct number of counter ions (Na + or Cl -) was added for neutralizing the total charge of the system followed by solvation of the system using TIP3P water molecules in a rectangular box with the buffering distances set to 12 Å in all directions. After that, the solvated system was submitted to an initial energy minimization process by applying Sander module (500 steps of steepest descent followed by 500 steps of conjugate gradient) followed by a 50 ps heating step where the temperature was gradually increased from 0 to 300 °K. After 50 ps of density equilibration, 500 ps of constant pressure equilibration at 300 °K with a time step of 2 fs was performed. Only bond lengths involving hydrogen atoms were constrained using the SHAKE algorithm. Final molecular dynamic simulation was individually performed for a range of 1 to 10 ns by applying the Particle Mesh Ewald (PME) method to calculate long-range electrostatic interactions. All calculations were done under periodic boundary condition where no constraint was applied to either the protein or the ligand molecules. The trajectory of the dynamic simulation was achieved by writing out the coordinates every 10 ps. After molecular dynamic simulation on receptor-ligand complex, snapshots were taken from the molecular dynamic trajectory with an interval of 10 ps. The dielectric constant values were set to 1.0 and 80 for the interior of solute and the surrounding solvent, respectively. Binding free energy was calculated for ligand-receptor complex using MM-GBSA. 80 The interaction energies for the snapshots were calculated while excluding water molecules and counter ions and presented as the average value in the RON-MSP system.

Results and Discussion
Binding free energies for the complexes of RON tyrosine kinase receptor and its ligand (i.e. MSP) as well as the mutants of either receptor or ligand (Table 1) were calculated by applying MM-GBSA method on molecular dynamic simulation data collected at different time. For this purpose, firstly the binding free energy was calculated for the RON-MSP wild type complex, then the residues involved in RON and MSP interaction were mutated to alanine followed by molecular dynamic simulation and re-calculation of the binding free energy for different time intervals ranging from 1 to 10 ns to estimate the contribution and the effect of individual residues in RON and MSP binding. The binding free energies (ΔG bind ) for wild type and mutant forms were calculated as follows: where G water (complex), G water (receptor), and G water (ligand) denote the free energies of the complex, receptor, and ligand, respectively. The free energy (ΔG) for each term is calculated using following equation:  Figure 1 shows the results of binding free energy calculations for the complex of wild type RON and MSP and their mutant forms using MM-GBSA method applied to molecular dynamic simulations ranging from 1 to 10 ns. These results have been also illustrated in Table 2. Results for ΔΔG binding (ΔG Binding-wild type -ΔG Binding-mutant ) for RON and MSPare also available in Table 3. The details of all calculations for mutants and wild types of receptor and ligand are available in appendices 1 and 2. Cancer is one of the most important causes of death in the world 83 and several strategies including pharmacotherapy protocols are employed to control this devastating condition. 84 Due to the importance of protein-protein interactions in cancer initiation and development, many efforts have been dedicated to target cancer cells by inhibition of those PPIs involved in cancer progression. 5,85,86 RON a tyrosine kinase receptor has gained considerable attention as promising target in cancer therapy. 68 Most of the therapeutic agents developed so far against RON interfere with RON and MSP binding highlighting the importance of PPIs. 69 Therefore, the identification of hot spots involved in the interface of RON-MSP complex is of great importance in rational drug design.
In the current study, the residues reported to be involved in RON-MSP interactions ( Figure 2) were virtually mutated to alanine one at the time to determine the contribution of each residue using MM-GBSA approach. The binding free energy difference between the mutant and the wild type complexes was obtained as follows:

ΔΔGBinding = ΔGBinding-wild type -ΔGBinding-mutant
In this expression, negative ΔΔG binding value implies that the substitution of the corresponding amino acid with alanine is an unfavorable substitution whereas a positive value indicates a favorable substitution in terms of binding free energy compared to the wild type complex. 87    The results of molecular dynamic simulation of RON indicated that all receptor (except for Glu 190 and Ser 195 ) and ligand (except for Arg 639 and Glu 644 ) mutants have low affinity compared to the wild type as deduced from the negative ΔΔG values shown in Figure 2 and Table 3.
One of the crucial residues at the interface of RON-MSP complex is RON Gln 193 ; its side chain NH 2 group makes two ionic interactions with carboxylate group of MSP Glu 644 and carbonyl group of Arg 639 . In addition, Arg 639 of MSP is involved in another interaction with RON Glu 190 which will be discussed later. 76 The MM-GBSA based binding energy calculations on the wild type and Q 193 A mutant showed that this amino acid is important in the binding (also confirmed by Chao and coworkers) 76 Table 3). It seems that the importance of these residues is related to their interactions with more than one residues on MSP (except for Pro 288 ). In the case of Pro 288 , it interacts only with MSP Arg 521 which in turn is highly important due to its participation in multiple interactions with RON Glu 287 , and Glu 289 . 76 Based on binding ΔΔG values, His 424 seems to be less important in comparison to other RON residues at the interface. However, this residue can also be considered as a hot spot on RON (Table 2 and Table  3). Such disagreement between the reported experimental results and our in silico estimates may be due to the fact that Glu 190 interacts with two different MSP residues (i.e. Arg 639 and Arg 683 ), which are already interacting with other RON residues. 76 Therefore, lack of their interactions with Glu 190 may not contribute favorably in the overall binding energy. RON Ser 195 is shown to be involved in a charge-charge interaction with MSP Arg 683 , 76 however, our results did not identify this amino acid as an important residue in RON-MSP complex (Table 3). Again the disagreement between our in silico estimates and the crystallographic data may be due to the formation of another interaction by Arg 683 with RON via Glu 190 which renders the interaction between Ser 195 and Arg 683 less important. 76 The only previous experimental site directed mutagenesis studies on the residues at the interface of RON-MSP complex was carried out for Arg 683 and results obtained are in agreement with the ones discussed below. 75 This amino acid is an important residue in the interaction of RON-MSPβ complex based on in silico calculation despite the results obtained for its partners on MSP (i.e Glu 190 and Ser 195 ). In order to gain more information regarding these residues an in silico double mutation (E 190 A/S 195 A) study was performed. This double mutation lead to a positive ΔΔG value indicative of their harmonic interplay in the interaction with Arg 683 .

Conclusion
In modern drug design and discovery process, computational approaches have streamlined a promising perspective by supplying useful and supportive information. In this context, identification of hot spots in biomolecules' interactions through the estimating the binding affinity of molecules towards targets of interest can provide valuable information where protein-protein interactions are important initiators in cancer pathogenesis. Virtual alanine scanning mutagenesis is one of the tools that are commonly employed for this purpose. Therefore, in the current investigation, amino acids reported to be at the interface of RON-MSP complex were evaluated using the MM-GBSA method and some of them were assigned as hot spots in the interaction. Taken