Effect of STAT3 inhibitor in chronic myeloid leukemia associated signaling pathway: a mathematical modeling, simulation and systems biology study

Chronic myeloid leukemia (CML) is a hematopoietic stem-cell disorder which proliferates due to abnormal growth of basophil cells. Several proangiogenic molecules have been reported to be associated in CML progression, including the hepatocyte growth factor (HGF). However, detail mechanism about the cellular distribution and function of HGF in CML is yet to be revealed. The proliferation of hematopoietic cells are regulated by some of the growth factors like interleukin 3 (IL-3), IL-6, erythropoietin, thrombopoietin, etc. In this study IL-6 pathways have been taken into consideration which induces JAK/STAT and MAPK pathways to decipher the CML progression stages. An attempt has been made to model these pathways with the help of ordinary differential equations (ODEs) and estimating unknown parameters through fminsearch optimization algorithm. Some of the specific component like STAT3, of the pathway has been analyzed in detail and their role in CML progression has been elucidated. The roles of STAT3 inhibitors into the treatment of CML have been thoroughly studied and optimum concentration of the inhibitors have been predicted. Electronic supplementary material The online version of this article (doi:10.1007/s13205-015-0357-7) contains supplementary material, which is available to authorized users.


Introduction
Chronic Myeloid Leukemia (CML) is a stem cell disorder resulted by abnormal growth of granulocytes (Strife and Clarkson 1988). The CML progression can be further differentiated into three phases namely Chronic, Accelerated and Blast depending on its severity, respectively (Kantarjian et al. 2006). The occurrence of CML has been reported to increase in two folds during last decade (An et al. 2010). It has been further reported that, translocation of the ABL gene on chromosome 9th to BCR gene on chromosome 22nd results into a hybrid Philadelphia chromosome with chimeric BCR-ABL genes (Daley et al. 1990;Rowley 1973;Ye et al. 2014). These chimeric genes encode BCR-ABL oncoproteins which contains the activated tyrosine kinase region of ABL, causing aberrant tyrosine kinase activity and hence promotes the proliferation of CML cells (National Comprehensive Cancer Network Clinical Practice Guidelines in Oncology, Chronic Myelogenous Leukemia, v2, 2009). The constitutive proliferations of BCR-ABL genes and their aberrant expressions lead to several abnormalities in the peripheral blood and bone marrow which consequently alters the number of granular leukocytes (Daley et al. 1990;Baccarani et al. 2006;Cerny-Reiterer et al. 2012). It is being reported that the earlier treatments processes involving chemotherapy, bone marrow and stem cell transplantation are often associated with serious side effects (Baccarani et al. 2006). Nowadays, CML is being treated by the BCR-ABL tyrosine kinase inhibitors, which leads to the increased rate of patient survivability. Most of inhibitors target the inhibition of ATP binding process with tyrosine kinase domain of ABL protein (Cerny-Reiterer et al. 2012;Zhang et al. 2010). The breakthrough treatment for CML was reported with the discovery of Imatinib (Savage and Antman 2002;Pray 2008). However, the reoccurrence of the CML cases in Imatinib treated patients has been reported due to the point mutation of amino acid residues in BCR-ABL receptors active site. Although few other receptors like STAT3, STAT5, and SRC have also been investigated as alternative to overcome the Imatinib resistance. In this study STAT3 receptor protein has been taken into consideration to explore the CML associated biochemical pathways like IL-6, which induces the MAPK and JAK-STAT (Reynaud et al. 2011;Bhalla and Iyengar 1999). IL-6 is an interleukin which is encoded by the IL-6 gene and responsible for pro-inflammatory cytokine and anti-inflammatory myokine activity. In overall IL-6, JAK/STAT and MAPK pathways, common gp130 (glycoprotein 130) binds to the plasma membrane receptor complex. Signal transduction activates the JAK tyrosine kinase family members, which activates the transcription factors of the STAT (Peifer et al. 2014). The interlinked pathway for IL-6-type cytokines is the MAPK pathway. Schematic representation of whole IL-6 associated signaling pathway is shown in Fig. 1. It has been reported that activation of STAT3 by bone marrow cells protects the CML from tyrosine kinase mediated inhibition of BCR-ABL protein. So, targeting the STAT3 in addition to BCR-ABL would enhance CML cells with kinase-independent resistance to TKI therapy. Recent study has reported BP-5-087 molecule as a multi-target inhibitor for BCR-ABL and STAT3 receptor. Further it was found that BP-5-087 molecule act on the Imatinib resistant cells and reduced the CML progression significantly (Eiring et al. 2015).

Methods
The entire experiments were performed on Centos 6.5 of Linux operating systems with 12 GB RAM, NVIDIA graphics and 1 TB Hard disk computer system. Hardware of the system was as follows: Intel(R) Core (TM) i7-3770 CPU @ 3.40 GHz processor.

Model description
CML associated signaling pathway in which IL-6 induces the two pathways, i.e. JAK/STAT and MAPK have been taken as model pathways. IL-6 is the first component of the said pathway, which enters into the cell by cell membrane and activates the entire pathway. In this process, IL-6 binds to the protein gp80, which is not involved in signaling process of the protein. The receptor gp130 which is a part of the signaling process attaches to the tyrosine kinases of JAK. Then IL-6-gp80 complex reacts to the compound gp130-JAK to form the complex IL6-gp80-gp130-JAK. Formation of this complex results in the dime formation of the complex IL6-gp80-gp130-JAK to form (IL6-gp80-gp130-JAK)2 (Heinrich et al. 2003). SHP-2 binds to the phosphorylated dimer complex (IL6-gp80-gp130-JAK)2 to form (IL6-gp80-gp130-JAK)*2-SHP-2 complex which act as the point of initiation for both the MAPK and JAK-STAT pathway (Yamada et al. 2003). In JAK-STAT signaling process, phosphorylated dimer (IL6-gp80-gp130-JAK)*2 recruits the transcription factor STAT3. The phosphorylation of the component STAT3 leads to its dissociation from the (IL6-gp80-gp130-JAK)2 complex and further undergo the process of dimerization. The dimerized complex of the STAT3 enters into the nucleus to undergo dephosphorylation and transferred back again into the cytosol for the another cycle of (Heinrich et al. 2003). It has been reported that apart from BCR-ABL, STAT3 receptor also plays major role in CML progression which is an important component in JAK/STAT pathway. We assumed the initial concentration of STAT3 in the pathway to be 1 nM as earlier reported (Charusanti et al. 2004). From the literature, it is been reported that BP-5-087 can inhibit the activity of STAT3 by 70 %, therefore we took the inhibited (treated) concentration of STAT3 as 0.3 nM (Eiring et al. 2015). JAK/STAT signal transduction pathway is regulated by the two phosphatases of STAT3 namely PP1 and PP2. Theses phosphatases enhance the activation of STAT3 component in nucleolus and cytosol both. The phosphatase PP2 present in the nucleus is a very important part of JAK-STAT since it is responsible in deactivating the phosphorylated dimer complex of STAT3 in the nucleus. This process of deactivation, play a crucial role for coming back of STAT3 to the cytosol so that it can undergo another cycle of phosphorylation-dephosphorylation. Binding of SHP-2 to the complex (IL6-gp80-gp130-JAK)*2 plays a major role in linking JAK-STAT to the MAPK pathway. This process of complex formation is an important process for the MAPK signaling pathway.

Kinetics of biochemical reactions
Schematic representation of chemical process of the pathway is shown in Fig. 1. The whole chemical processes in the pathway are decomposed into elementary reactions. This reaction model contains 67 state variables and 111 unknown rate constants. It has been assumed that initial concentration of each component is 1 nM and experimental concentrations are taken from reported literature (Charusanti et al. 2004). The kinetic reactions of the pathway are provided in Appendix I (Online Supplementary Material).

Mathematical modeling of JAK/STAT and MAPK pathway
In this work, model parameters were calculated from the experimental data obtained from IL-6 induced signaling pathways. Systems biology techniques were used to understand various components involved in reaction mechanism of this pathway and to predict the behavior of individual reaction components (Heinrich et al. 2002;Papin et al. 2005;Tyson et al. 2003;Ge et al. 2013). Deterministic rate laws have been used to simulate various chemical and biochemical reactions (Edda et al. 2006;Tian and Song 2012;Klipp and Liebermeister 2006). The deterministic mathematical modeling of chemical network works on the basis of law of mass action and Michaelis-Menten reaction kinetics (Kanodia and Finn 2014). It has been further reported that the prior information of initial concentration, law of mass action etc., can enable the simulation of time series data to study substrate concentration and its temporal behavior (Edda et al. 2006;Tian and Song 2012). However, simulation process provides in depth understanding of each components of the biological systems and its activities in different environmental conditions (Gilbert et al. 2006;Guerriero et al. 2009). The rate of change in substrate concentration over a given period of time can be utilized to find the crucial parameters involved in the reaction mechanism (Salis and Kaznessis 2005). IL-6 induced signaling pathway has been simulated through deterministic mathematical modeling (Kwiatkowska et al. 2006). The candidate reaction pathways were simulated through ordinary differentiation equations to get the individual rate constant parameters (Deng and Tian 2014;Liepe et al. 2014). The parametric optimization was performed by ''fminsearch'' algorithm on MATLAB R2012a. Complete methodology has been shown in Fig. 2 In this process IL-6 induces two pathways like JAK/ STAT and MAPK pathway which are responsible for various cellular processes.

Parameter estimation
The parameters associated with above stated JAK/STAT and MAPK pathway's reactions are estimated through ODE solver. This numerical simulation of model was performed by using MATLAB ODE solver toolbox (Bogacki and Shampine 1989). The standard ODE solver ODE45 was used to construct and solve the differential equations based on Runge-Kutta approximation method. Runge-Kutta methods involves implicit and explicit iterative methods for temporal discretization and approximation of ODE solutions. The proposed models with continuous states were being iterated and further optimized by fminsearch.
where '@function' is the vector of function handles, 'tspan' a vector specifying the time span of simulation and 'x0' is the state variables of differential equations. The 'tspan' and 'x0' represents the specific time interval and initial starting conditions, respectively.

Parameter optimization
fminsearch is a MATLAB function which minimize the error function with multiple variables, with the assigned estimated values of state variables. It utilizes nonlinear unconstrained multivariable optimization problems based on Nelder-Mead Sequential Simplex (NMSS) algorithm (Andrews and Larry 1992; Lagarias et al. 1998).
With the multiple functions are f(x1), f(x2)… f(xn ? 1) the Nelder Mead Algorithm. NMSS is a nonlinear optimization technique, which is used to solve numerical problems with unknown derivatives. However, NMSS is a heuristic search method that can converge to non-stationary points (Powell 1973).
where function is an M-file function such as whereas, @minimize is a function containing objective function. Error function used for minimizing scalar function is: Experimental data Experimental concentrations of the each substrate have been taken from (Singh et al. 2006) study. All the experimental values are shown in Appendix 1 (Online Supplementary Material) which is required for the modeling and simulation.

Results and discussion
It has been reported that CML exhibit an aberrant expression of beta lymphoid cells which are directly influenced with cytokine IL-6 pathway, which involves two significant pathways like JAK/STAT and MAPK pathway. To get insight of the molecular mechanism of the CML, extensive in silico study of both pathways are needed. Aberrant expression of STAT3 in the IL-6 signaling pathway had been proven as prime cause of Chronic Myeloid Leukemia (CML). It has been established that STAT3 is a potential therapeutic target for the treatment of CML because it is an important component of the JAK/STAT pathway. Alteration or inhibition of STAT3 concentration indirectly influence the another associated MAPK pathway.

Mathematical model of the JAK/STAT and MAPK pathway
These two interlinked pathways consist of 67 biochemical reactions and their 111 unknown rate constants. These biochemical reactions involve mass-action kinetic and Michaelis-menten kinetic reactions. All these reactions of the MAPK and JAK/STAT pathways were mathematically modeled in the form of ordinary differential equations. In continuation, all differential equations were solved by ode45 solver in MATLAB. ODE models of the said chemical reactions are given in Supplementary material 3.
In this work role of STAT3 has been analyzed, as STAT3 is the most significant component of JAK/STAT pathway. STAT3 combines with (IL6-gp80-gp130-JAK) complex and form another complex, i.e. (IL6-gp80-gp130-JAK)-STAT3. STAT3 activates the (IL6-gp80-gp130-JAK) complex to enter into the MAPK pathway by dissociating themselves before the entry. By analyzing complete pathway of JAK/STAT and MAPK one can infer that STAT3 plays a decisive role into CML associated pathway.

Parameter estimation
Experimental concentrations of the substrates involved in IL-6 induced signaling pathway are shown in Appendix I  Table 1.

Model simulation
The mathematical model of 67 state variables has been simulated with respect to time for a period of 32 min to understand the behavioral changes of the various components involved in the signaling pathway. The simulation results consists of both, the normal condition and the drug treated condition for all the components. The significant changes in the components have been reported in the following simulated results.

STAT3 complexes
Existence of Cytosolic STAT3 is important for JAK/STAT signaling pathway which mainly involves the transcription factor, STAT3 for binding with various complexes in the pathway. If STAT3 is inhibited, it can hamper the signaling pathway through JAK/STAT and MAPK pathway. The Figs. 3a, b, and 4a, b represents the behavior of STAT3 complexes in normal condition and treated condition. Inhibition of STAT3 leads to decrease in the concentration of STAT3 complexes, such as (IL6-gp80-gp130-JAK)*2 -STAT3C, STAT3C*, (IL6-gp80-gp130-JAK)*2- STAT3C* and STAT3C*-STAT3C* as compared to the normal condition. Phosphorylation of STAT3 is essential for the signaling through JAK-STAT. These complexes help in the phosphorylation and transportation to the nucleus. If the STAT3 is inhibited, then it cannot be phosphorylated and consequently the JAK/STAT pathway will be inhibited.

PP1 and PP2 complexes
JAK-STAT signal transduction pathway is regulated by the two STAT3 phosphatases, PP1 and PP2. Theses phosphates influence the activation of STAT component in both the cytosol and the nucleus. Figures 6b and 7a represents the normal and treated condition of the PP1 complex, PP1-STAT3C* and PP1-STAT3C*-STAT3C*. Inhibition of STAT3 also influences. Inhibition of STAT3 component affects the PP1 and due to which the concentration of PP1 complexes are also reduced. If PP1 is affected, it will hamper the activation of STAT component then further the JAK/STAT pathway. Along with STAT3, PP2 also plays a crucial role in JAK/STAT signaling pathway. Figure 7b represents the behavior of PP2 with respect to time when the STAT3 is in normal and treated condition. It can be analyzed that level of PP2 in treated condition is reaching to the lower level as compared the normal condition at the end of the simulation. This shows that STAT3 is affecting the levels of PP2 in the pathway. PP2 deactivates the phosphorylated STAT3 dimer STAT3 N-STAT3 N* in the nucleus. This deactivation of STAT3 is necessary for the phosphorylation and dephosphorylation cycle of STAT3 in the JAK/STAT signaling pathway. If the STAT3 is inhibited, the level of PP2 will be decreased and it will not be able to deactivate the phosphorylated STAT3 dimer in the nucleus. This condition will affect the concentration of STAT3 in cytoplasm and STAT3 N*-STAT3 N* in the nucleus and will further reduce the efficiency of signaling process through JAK-STAT pathway and other cellular processes.

Conclusion
In this work, mathematical models of all biochemical reactions of CML associated signaling pathways like JAK/ STAT and MAPK have been generated. An attempt has been made to investigate the effect of inhibitors on particular receptor like STAT3 receptor in CML associated signaling pathways. Each step of pathways was converted into respective biochemical kinetics following mass action principle of michaelis-menten kinetic reactions. We proposed an effective optimization approach for estimating model parameters of the JAK/STAT and MAPK signaling pathways. Fminsearch function was used to minimize the error during simulation process. It generates observation data at different time interval together with the first order derivatives of the model parameters. The calculated parameters values were used to infer unknown parameters. Simulation results suggested that the proposed optimization method is effective, robust and reliable parameter estimation technique. Further the role of STAT3 has also been predicted in CML associated signaling pathways. Since the phosphorylation of STAT3 is crucial for signaling process of JAK/STAT and MAPK pathway, hence its inhibition affects other components of the pathways such as PP1, PP2, and SHP2. Inhibitory effect of the BP-5-087 ligand on STAT3 receptor was compared with the normal (untreated condition). It can be inferred from the simulation results that BP-5-087 can successfully reduce the progression of CML which is further in concurrence with the earlier reported in vitro study by Eiring et al. (2015). The proposed model may be utilized to analyze an overall impact of various inhibitors targeting the CML associated receptor proteins.