Parameters Involved in Autophosphorylation in Chronic Myeloid Leukemia: a Systems Biology Approach

Background: Chronic myeloid leukemia (CML) is a stem cell disorder characterized by the fusion of two oncogenes namely BCR and ABL with their aberrant expression. Autophosphorylation of BCR-ABL oncogenes results in proliferation of CML. The study deals with estimation of rate constant involved in each step of the cellular autophosphorylation process, which are consequently playing important roles in the proliferation of cancerous cells. Materials and Methods: A mathematical model was proposed for autophosphorylation of BCR-ABL oncogenes utilizing ordinary differential equations to enumerate the rate of change of each responsible system component. The major difficulty to model this process is the lack of experimental data, which are needed to estimate unknown model parameters. Initial concentration data of each substrate and product for BCR-ABL systems were collected from the reported literature. All parameters were optimized through time interval simulation using the fminsearch algorithm. Results: The rate of change versus time was estimated to indicate the role of each state variable that are crucial for the systems. The time wise change in concentration of substrate shows the convergence of each parameter in autophosphorylation process. Conclusions: The role of each constituent parameter and their relative time dependent variations in autophosphorylation process could be inferred.


Introduction
Chronic myeloid leukemia (CML) is a hematopoietic stem cells disorder caused by translocation of 9th chromosome of ABL oncogene to 22nd chromosome of the adjacent BCR oncogene (Bayard et al., 1988;Daley et al., 1990). The specific chromosome where translocation occurs is known as Philadelphia chromosome which is present throughout the CML diseased cells (Eaves C et al., 1993). Aberrant expression of the BCR-ABL oncogenes resulted in to BCR-ABL oncoprotein which is found in more than 98% of the CML cases Mauro et al., 2001). It has been reported that Imatinib acts as an efficient inhibitor to prevent the fusion of BCR-ABL oncoprotein under the target based therapy (Marley et al., 2000). This protein acts as a kinase enzyme which aberrantly enhances the addition of phosphate group to the SH2 domain of BCR-ABL protein substrate. It has been further reported that the activity of the BCR-ABL kinase protein is dependent upon extracellular and intracellular signaling processes (Clarkson et al., 1991, Xi-Shan et al., 2014. The expression of these enzymes regulates the on/off mechanism of the cellular proliferation process in CML (Clarkson et al., 1991). Several studies focusing on the CML pathway especially targeting the tyrosine kinase inhibitors like Imatinib, Baustinib, Dasatenib etc.,were being reported earlier (Clarkson et al., 1991;Yuan-Xin et al., 2014). In this study, overall mechanism of autophosphorylation of BCR-ABL was analyzed through mathematical modeling and simulation to understand the roles of individual component in proliferation of CML.
Autophosphorylation is the post translational modifications process where the phosphorylation of the enzyme kinase is resulted by addition of phosphate group to specific amino acid residues (serine, threonine and tyrosine etc.), to activate the catalytic activity (Summers, 2011). This process can be further classified into cis-autophophorylation (when a kinase's own active site catalyzes the phosphorylation reaction) and transautophosphorylation (when another kinase of the same type provides the active site that carries out the reaction often during dimerization) (Smith et al., 1993).
The onset of CML disease hampers hematopoietic stem cells hence resulting in to aberrant hematopoiesis i.e. abnormal blood formation (Kagita et al., 2010). The detailed study of hematopoiesis is of substantial importance to understand the role of autophosphorylation process of the CML. It has been reported that BCR-ABL has a unique ATP binding site SH-2 domain near to the substrate proteins binding region (Savage et al., 2002).The details about BCR-ABL autophosphorylation pathway, its associated receptors and their binding patterns are shown in Figure 1.
Various systems biology approaches have been reported which utilized the mathematical modeling techniques for identifying the regulatory mechanism of biochemical pathways (Obel et al., 2014;MacLean et al., 2015;Pokhilko et al., 2015;Vera et al., 2015). These mathematical models were successfully implemented to analyze various interactions among components of the particular pathways and simulate them to the environmental factors or external signals (Marshall-Colon et al., 2014;Finley et al., 2014;Wu et al., 2014). On the other hand, it also enables to explore the system function at various levels and help to generate hypothesis about biological experiments. However, there are some limitations of this method, like identification of unknown parameters which are estimated from experimental results (Liepe et al., 2014;Denget al., 2014;Raue et al., 2014).

Materials and Methods
In this study, model parameters were calculated from experimental data obtained from autophosphorylation of Bcr-Abl pathway of chronic myeloid leukemia cells. Systems biology techniques were used to understand various components involved in reaction mechanism of Bcr-Abl pathway and to predict the behavior of individual reaction components (Heinrich et al., 2002;Papin et al., 2005;Tyson et al., 2003;Pep et al., 2004;Yu-Qing Ge at al., 2014).It utilizes the deterministic rate laws to simulate various chemical and biochemical reactions (Edda et al., 2006;Tianhai et al., 2012). The deterministic mathematical modeling of chemical network works on the basis of law of mass action and Michaelis Menten reaction kinetics. 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;Tianhai et al., 2012). However, simulation process provides in depth understanding of each components of the biological systems and its activities in different environmental conditions. 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 (Howard et al., 2005). Autophosphorylation reaction mechanism of the CML has been simulated through deterministic mathematical modeling. The candidate reaction pathways were simulated through ordinary differentiation equations to get the individual rate constant parameters. The parametric optimization was performed by "fminsearch" algorithm on MATLAB R2012a.
The entire simulation was performed on Centos 6.5 of Linux operating systems with 12 GB RAM, NVIDIA 2GB graphics card and Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz processor.
Our model contains autophosphorylation process of Bcr-Abl and representation of this process is constructed in Figure 2. In this process, unphosphorylated Bcr-Abl substrate binds to ATP to form a Bcr-Abl·ATP complex. This Bcr-Abl·ATP complex then binds with another Bcr-Abl·ATP complex. Complete complex formation leads to cross-phosphorylate Bcr-Abl and forms phosphorylated form of Bcr-Abl. This activated form of Bcr-Abl is represented as Bcr-Abl*.
The conventional modeling approach for signaling pathways involves solving individual chemical reactions through ordinary differential equations (ODEs). Whereas, individual chemical reaction can be modeled by decomposing the signaling pathway into elementary reactions. Law of mass action kinetics can be applied to each elementary reaction to obtain ODE's as shown in Eq.1. Reactions in this pathway are represented in Eq.12 to 15 using law of as mass-action kinetics. Hence the mathematical representations of each mass action kinetic reactions involved in autophosphorylation are as follows:

3) d[L]/dt= -ka[R][L]+ kd[RL] (4) d[RL]/dt= ka[R][L]-kd[RL] (5)
The reason for using mass-action kinetics in autophosphorylation signalling pathway is because of the value of substrate concentration is much larger as compared to the enzyme concentration. Mathematical modelling of mass action kinetic reaction networks using ordinary differential equations produced polynomial vector fields. Concentration of each component of the reaction is balanced throughout the pathway modelling. If a component is formed at any time of reaction then it is consumed in some reaction so as to maintain the overall equilibrium of substrate concentration in the pathway. In general, differential equation of a particular component (A) is written as:

Rate of change of [A] = Amount of [A] formed -Amount of [A] consumed (6)
Numerical simulation of model was performed by using MATLAB ODE solver tool box (Bogacki et al., 1989). The standard ordinary differential equations (ODEs) solver i.e. ODE45 was used to construct and solve the differential equations based on Runge-Kutta approximation method. Runge-Kutta methods are an important family of implicit and explicit iterative methods, which are used in temporal discretization for the approximation of solutions of ordinary differential equations. Therefore, ODE45 was selected as the default solver for the proposed models with continuous states.

[t,x] = ode45(@function,tspan,x0) (7)
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.
fminsearch starts at the point x0 and returns a value x that is a local minimizer of the function described in @ function. Error function used for minimizing scalar function is:

Results
To address the issue of insufficient experimental data of autophosphorylation process of Bcr-Abl pathway, we propose parameter estimation to make reliable inference of model parameters. Mass-action kinetic reactions involved in autophosphorylation process of Bcr-Abl, were mathematically modeled in the form of ordinary differential equations. Complete chemical process of autophosphorylation of Bcr-Abl is shown in Figure1, and decomposed themto elementary reactions. This reaction model contains four state variables and five unknown rate constants. It has been assumed that initial concentration of each component is 1 nM and experimental concentrations are taken from reported literature (Pep et al., 2004).
ODE models of the chemical processes are given below:
Experimental concentrations of the substrates involved in autophosphorylation process of Bcr-Abl, ATP, Bcr-Abl·ATP and Bcr-Abl·ATP•Bcr-Abl·ATP are 0.2, 7, 0.12 and 0.24 nM respectively. Where ode45 solves the differential equation models given in Eq-12 to Eq-15 and fminsearch minimizes the error function and finds the optimized values of rate constants kai and kdi. Optimization technique is use to get the as closure as possible to the reported experimental value. We initiated the simulation by considering unit value i.e. 1 of each component of the reaction for a period of 3 seconds. Since our model consists of four state variables, we divided simulation time into four intervals. After selection of random time intervals, 0 to 3 seconds give the most closure value to the reported one. The results of simulation and values of unknown optimized rate constants k a1 , k d1 , k a2 , k d2 , kautoare 0.1831, 1.3898, 0.8557, 0.7021, and 0.0033 respectively.
Plot analysis of components is as follows, [Bcr-Abl] concentration initially increases exponentially from 0s to 0.7s. After 0.7s the [Bcr-Abl] concentration gradually decreases as shown in Figure 3A. Concentration of [Bcr-Abl.ATP] initially decreases and then gets stable after 1second as shown in Figure 3B.

Discussion
This work presents an optimization approach to address the issue of inadequate experimental data for inferring unknown parameters in mathematical models of Autophosphorylation reactions. In this deterministic simulation process, mathematical modeling of chemical reactions involved in autophosphorylation of Bcr-Abl was done. Chemical reactions of the process were represented as mass-action kinetics by using ordinary differential equations (ODEs). The ODE45 solver was used to construct and solve first order coupled nonlinear differential equations for deterministic simulation and fminsearch to find the optimized value of unknown rate constants involved in these reactions. Mathematical modeling has been extensively used earlier for many biochemical reactions to make appropriate assumptions of rate constants. This is majorly useful for experimentally untested signaling pathways and helps to predict unknown rate constants. One of the major steps in developing mathematical models is to estimate unknown parameters  DOI:http://dx.doi.org/10.7314/APJCP.2015.16.13.5273 Parameters Involved in Autophosphorylation in Chronic Myeloid Leukemia: a Systems Biology Approach of the model based on experimentally measured quantities. Substrate, enzyme and product behavior can also be determined with respect to time. In this model of autophosphorylation of Bcr-Abl was done protein, substrate concentration of [Bcr-Abl·ATP•Bcr-Abl·ATP] increase with respect to time till 1 second, after that it starts plays hardly any role in the pathway. Simulation of [Bcr-Abl·ATP•Bcr-Abl·ATP] will not affect the autophosphorylation process of Bcr-Abl protein. All unknown parameter at which optimization techniques gives minimum error were calculated as k 1 (0.9097), k -1 (1.1023), k 2 (0.5196), and k 3 (-0.0448).
In conclusion, this work proposes an effective optimization approach for estimating model parameters in mathematical models of autophosphorylation process. This method uses an fminsearch function to approximate the underlying solution of the autophosphorylation model. The fminsearch function generates observation data at different time interval together with the first order derivatives of the model reaction. All generated information is used to infer unknown parameters from the reactions involved in autophosphorylation process. An attempt has been made to address the issue of inadequate experimental data for estimating unknown parameters. Simulation results suggest that the proposed optimization method is effective and robust approach reliable parameter estimation. It has been further concluded that interaction between Bcr-Abl oncoprotein with ATP molecules is crucial for the proliferation of cancerous growth. Hence it is inferred that, Bcr-Abl could be a putative target for CML prevention.