On Nonlinear Dynamics of the STAT5a Signaling Protein

In this paper we model dynamics of cross talk between MEK/ERK and JAK/STAT signaling pathways by means of nonlinear ordinary differential equations. The considered system of four ordinary differential equations is reduced to one ordinary differential equation, representing dynamics of the phosphorylated STAT5a signaling protein. We show that the diffusion together with the corresponding biochemical reactions is likely to play a critical role in governing the dynamical behavior of the considered signaling protein. By the modified method of simplest equation to the described reaction-diffusion equation we obtain an analytical solution which explains drop and jump propagation of the STAT5a protein concentration. Keywords-STAT5a signaling protein; PDEs; Modified method of simplest equation; analytical solution; drop and jump propagation;


I. Introduction
The idea to modeling proteins interactions in the intracellular environment by spatial-temporal systems is based on the following circumstance: While the dynamics of the interacting proteins and their molecular pathways and networks can be described by reaction models, the heterogeneous distributions of the protein concentrations are not taken into consideration.It is proved, however, that similar inhomogeneous protein distributions in the form of cellular jump or drop propagation play an important role in the control of the main processes in the cell.In this way the cellular complexity appears to be space-temporal, expressed mathematically by reaction-diffusion models.In fact, the traditional approximation scheme of a well-stirred reactor is most used and somewhat plausible simplification.However, the concentration gradients of cell enzymes that modulate signal transduction refute this simplification [1][2][3][4].The role of diffusion in reaction-diffusion systems of the cell becomes significant when reactions are relatively faster (but not so very) than diffusion rates.Sometimes the term 'molecular crowding' is used to denote more specific type of spatial distribution [5,6].The biochemical essence of this phenomenon is based on the fact that the state of phosphorylation of target molecules with spatially separated membrane-localized protein kinases and cytosolic phosphatases depends essentially on diffusion.Then the very high protein density in the intracellular space, or so called molecular crowding, mentioned above can enhance the spatial effect.Consequently, molecular crowd-ing can also alter protein activities and break down classical reaction kinetics [5].Therefore here we introduce reaction-diffusion modeling and its computational tools to concrete example of ERK (Extracellular signal-Regulated Kinase) and STAT (Signal Transducers and Activators of Transcription) protein interaction.In [7] the cross talk between the MEK/ERK and JAK/STAT signaling pathways has been modeled in a form of a system of nonlinear ordinary differential equations for the protein concentrations (JAK is an abbreviation of Janus Kinase, and MEK is an abbreviation of MAPK/ERK Kinase).Next, in the same paper the spatial modeling of the mentioned interaction has been performed, by introducing an appropriate diffusion-reaction scheme.By the accomplished stability analysis in [7] the authors concluded that in terms of the Turing bifurcation in ERK and STAT dynamical model the mentioned above crowding effect can be interpreted as a process of stabilization of the dissipative structures inherent to the considered intracellular process.

II. The method of the simplest equation
There are many approaches for obtaining analytic solutions of nonlinear partial differential equations [8][9][10].In this paper we will use the modified method of the simplest equation.The method is a modified version of the method of the simplest equation, created by N. A. Kudryashov [11,12] that is based on the fact that after application of an appropriate ansatz a large class of nonlinear PDEs can be reduced to nonlinear ODEs of the kind (P means polynomial) and for some equations of the kind (1) particular solutions can be obtained which are finite series constructed by the solution Φ(ξ) of a simpler equation referred to as the simplest equation.The simplest equation can be the equation of Bernoulli, equation of Riccati, etc.The substitution of (2) in (1) leads to the polynomial equation where the coefficients σ i , i = 0, 1, ..., r depend on the parameters of the equation and on the parameters of the solutions.Equating all these coefficients to zero, i.e., by setting one obtains a system of nonlinear algebraic equations.Each solution of this system leads to a solution of kind ( 2) of (1).In order to obtain a non-trivial solution by the above method we have to ensure that σ r contains at least two terms.To do this within the scope of the modified method of the simplest equation we have to balance the highest powers of Φ that are obtained from the different terms of the solved equation of kind (1).As a result of this we obtain an additional equation between some of the parameters of the equation and the solution.This equation is called a balance equation [13,14].We note that the method of the simplest equation and its modified version are closely connected to the problem for obtaining meromorphic solutions of nonlinear partial differential equations [15,16].By the methodology described in [15,16] one can obtain other interesting classes of solutions of nonlinear PDEs such as rational solutions for example.In addition we stress that by means of the traveling wave ansatz one reduces the nonlinear PDE to a nonlinear ODE and after this if an appropriate simplest ODE exists then a particular solution can be obtained that usually depends on as many parameters of the problem as possible.In many cases such particular solutions are among the few possible exact analytic solutions of the studied nonlinear PDE.

III. The interaction between ERK and STAT signaling pathways A. A reaction model and its reduction
The specificity of biological responses is often achieved in a combinatorial fashion through the concerted interaction of signaling pathways [17].Here we will discuss the eventual interaction between the units of MEK/ERK and JAK/STAT signaling pathways.In cell signaling the pathways are understood as networks of recurrent biochemical reactions, by which the information transmission (in a signal form) is accomplished.The disturbance of the intracellular signal transmission from the membrane receptors to the nucleus genes is assumed as a general reason for a cancer diseases progress.In [17] the authors discussed the basic features of interaction domains, and suggested that rather simple binary interactions can be used in sophisticated ways to generate complex cellular responses.Moreover, in [18], the protein STAT was found to play important roles in numerous cellular processes including immune responses, cell growth and differentiation, cell survival and apoptosis, and oncogenesis.The STAT target genes include SOCS/CIS, a class of inhibitory proteins that interfere with STAT signaling through several mechanisms.(SOCS is an abbreviation of Suppressor Of Cytokine Signaling and CIS means Cytokine Inducible SH2 domain containing).The protein SOCS/CIS can inhibit the STAT phosphorylation and block the access of STAT to receptors of JAKs or both [19].On the other hand, SOCS-3 can bind to sequester the Ras-GAP proten (The Ras proteins are members of the MEK/ERKpathway) [20].However, this is not only the way of interaction between STAT and ERK pathways.In [21][22], the authors suggested that the STAT5 functional capacity of binding DNA could be influenced by the Mitogen-Activated Protein Kinase (MAPK)-pathway.Later on, in [23] the interactions between STAT5a and ERK1 (or ERK2) signaling proteins was considered.A simple biochemical diagram of this interaction is presented in Fig. 1.According to Fig. 1 in [7] the following system of ordinary differential equations for the kinetics of STAT5a/S phosphorylation and ERK activation was constructed: Here e 1 , e 2 , s 1 and s 2 are state variables, representing concentrations of ERK-inactive, ERK-active, STAT5a-non-phosphorylated and STAT5a-phosphorylated proteins respectively.The following initial conditions are put on the system (5): They correspond to the initial concentrations of the above-mentioned proteins in mM units [27].Moreover, k 1 is proportional to the frequency of collisions of ERK and STAT5a protein molecules and presents a rate constant of reactions of associations; k 2 and k 3 are constants of exponential growths and disintegration; I is an inhibitor source respectively.The source I inhibits the phosphorylation of non-phosphorylated STAT5a.More concrete interpretation of the inhibitor I can be given by considering the role of thee SOCS proteins in linking the JAK/STAT pathway.Biological responses elicited by the JAK/STAT pathway are modulated by inhibition of JAK (and respective attenuation of STAT) by a member of the SOCS proteins [24,25].In [7] the inhibitor I is presented mathematically in the following manner: where Σ is a constant concentration of SOCS proteins and k 0 is a reaction rate constant of inhibition respectively.It is clear that the considered interaction between ERK and STAT pathways can occur only if Σ is sufficiently small, i.e. phosphorylation of the protein STAT5a exists.Therefore we assume the concentration of SOCS proteins to be sufficiently small.In addition in [7] the authors reduced (5) to the following two-dimensional system: 7) taking into account that only two equations of the four ones are independent.This means that between the concentrations e 1 .e 2 , s 1 and s 2 the following relations exist: where E 0 and S 0 are initial values of the sums of corresponding concentrations of inactive and active ERKs and non-phosphorylated and phosphorylated STATs.We assume for (7) that the inequality e 2 << s 2 holds, in view of the fact that the amount of ERK molecules is essentially smaller than the amount of STAT5a ones [21][22][23].In this way the inequality e 2 /s 2 = << 1 will be valid, and ( 7) can be written as: Next, in accordance with the QSSA (Quasi-Steady-State Approximation) theorem [26] we consider the first equation of ( 9) to be linear with respect to e 2 and we treat it as an attached system [26], i.e. we take into consideration that e 2 is a sufficiently small "constant".Further according to the requirements of the QSSA theorem we prove that the attached system has a stable steady state (then well-known Lyapunov definition of stability is satisfied).After replacing the steady state value in the second equation of ( 9) (the degenerate system), the quasi-stationary approximation of (7) (or ( 5)) is obtained in the form Next, for our convenience we rewrite (11) in the following manner: where we denote the variable s 2 only by s.
According to (11) the new coefficients are positive because we assumed that the initial concentration of the STAT5a protein is large, the ERK concentration is smaller than the STAT concentration, but bigger than the concentration of the SOCS proteins, and k 0 , k 1 , k 2 and k 3 are positive constants.

B. Stability analysis of the steady state solution of (12)
The fixed points of ( 12) can be found by solving the nonlinear equation They correspond to the following stationary concentrations s 0 : at Ds − E 0, which are positive, because and if the SOCS concentration Σ is sufficiently smaller than the total ERK (E 0 ) and STAT (S 0 ) concentrations as we assumed in the previous paragraph.According to the nonlinear character of the first term of ( 14) we can present it in the following form: The new coefficients C 1 and C 2 can be found by the equality or Moreover, the following relationship among the coefficients of the basic form exists too.
We replace the coefficients (20) in (18) and put it in (14).In this way we obtain the following steady state of ( 12): which is positive in view of ( 16) and ( 17).However, it will be valid only if (21) is satisfied.Taking into account (13), the last condition will hold if the coefficients k 1 or k 2 are sufficiently small (almost approaching zero), i.e. associations between ERK and STAT5a protein molecules or inactivation of active ERK molecules are very slow processes.
Let us now analyze the stability of this equilibrium.For the purpose we postulate the substitution s = s 0 + η, where η is a small variation (perturbation) from the steady state value s 0 The corresponding variational equation has the form: where It is easy to show that for s 0 = s 0 1,2 the coefficient w will have negative value.This is a sufficient condition for verification of the asymptotically stable character of the corresponding equilibrium states.On contrary, for s 0 = s 0 3 , w will have positive value, i.e. the corresponding equilibrium state will be unstable.Therefore we proved that the smallest (s 0 1 ) and biggest (s 0 2 ) steady states are stable, and the intermediate steady state s 0 3 is unstable, respectively.
Finally, the results in this section can be summarized as follows: 1) The dynamical behavior of the interaction between ERK and STAT5a proteins (in particular, between MEK/ERK and JAK/STAT signaling pathways) near its quasi-stationary state determines only by the behavior of the STAT5a signaling protein; 2) Stabilization of the considered process in time will observe at smallest (about 10 −3 mM ) and biggest (over 10 −2 mM) STAT5a concentrations.Destabilization of the same process will occur at intermediate steady state concentrations of the STAT5a protein, which will lead to initiation of some pathological process (for example, cancerogenesis).

IV. Space-temporal modeling of the STAT5a signaling protein A. A reaction-diffusion model of the STAT5a signaling protein. Stability analysis of the inhomogeneous distribution of the STAT5a concentration.
We introduce systems with distributed variables (reaction-diffusion models) when the connections between neighbor points of the space are taken into account by the diffusion law of molecular motion from the higher to lower concentrations.Here we take into account the reaction-diffusion effect, described in the introduction of this article to the reaction model of the STAT5a signaling protein.As a result we obtain the following one-dimensional model with distributed parameters where r is the space coordinate from the cell membrane to the nucleus, and D s is the diffusion coefficient of the STAT5a concentration.In order to analyze qualitatively and solve quantitatively (25), it is necessary to fix initial and boundary conditions for the function s(r, t) in the form: where l is the distance between the membrane and nucleus.Next, we consider the equation ( 25) and search for solutions of the kind of traveling waves: s = s(ξ) = s(r − vt) , where v is the velocity of the wave.In this way, by introducing the travelingwave coordinate ξ = r − vt, the equation ( 25) transforms to the following ordinary differential equation of second order: and next -in the following system of ordinary differential equations of first order: where denotes d/dξ.By analogy with the previous paragraph the fixed points of ( 27) can be found by solving the equation ( 14).Next we will analyze stability of the fixed points of (27) (in particular of (26) or (25)) in the phase plane (s, s').For the purpose we introduce the substitutions s = s 0 + and y = y 0 + ζ, where and ζ are small variations (perturbations) from the inhomogeneous equilibrium (s 0 , s 0 ).We substitute the last expressions in (27), and obtain the following variational equations: where w is a coefficient, presented by (24).The corresponding characteristic equation has the form: where are Routh-Hurwitz coefficients.According to the RouthHurwitz stability criterion they must have positive values to ensure the stable character of the corresponding steady state.Here the diffusion coefficient D s must be positive and the wave velocity v must be different from 0. Thereby the stability of inhomogeneous steady state of ( 27) (or ( 26)) is determined by signs in front of the coefficients v and w.We can observe a stable steady state solution of ( 27) if v > 0. If v < 0 the corresponding steady state will be always unstable.Thus, at v > 0 the smallest (s 0 1 ) and biggest (s 0 2 ) steady states are unstable in view of the negative values of w, and the intermediate steady state s 0 3 is stable because w is positive.
Finally, we can conclude that the average steady state concentrations of the STAT5a protein are indicative for the dissipative structures existence, but too low and too high ones are not.In this sense, the well-known crowding effect will hold, and only increase or decrease of the steady state distributions could assure the dissipative structure stability in the cell.

B. Traveling wave solutions of (26). Application of the modified method of simplest equation
The equation (26) (in particular (25)) will have solutions for such values of s, at which its nonlinearity can be reduced to a polynomial nonlinearity.The reduction of a non-linearity to a weak non-linearity is a commonly used approach in the nonlinear dynamics.In this way, here we develop the last three terms of (26) in a Taylor series centered in s = 0. We retain only the terms up to cubic power and obtain the following approximation of (26): where the new coefficients have the form: In this case we have not interested in the concrete values of the above-given coefficients, because their signs (positive or negative) do not affect our further analytical investigations.In addition, the cubic polynomial approximation means that we accept a weak non-linearity (but not linearization) of ( 26), i.e.D (the rate constant k 1 ) is sufficiently smaller than E (the rate constant k 2 ) to assure the approximation validity.Indeed, the last inequality follows from the biochemical consideration that the processes of ERK inactivation and STAT dephosphorylation are faster than that of ERK and STAT interaction.The last is of molecular recognition type [21][22][23].Next we will consider the wave propagation of the STAT5a density.For the purpose we will apply the methodology from Section II to (31).In order to solve (31) we constrict a solution as finite series where φ is the solution of a simpler equation (referred to as simplest equation according to Section II), and a i and c j are parameters that we will determine below.We substitute (33) in (31) and obtain an equation that contains powers of φ.
Next, we balance the highest power arising from the second derivative in (31) with the highest power arising in the term containing s 3 in the same equation.The resulting balance equation is r = 2n + 2, n = 2, 3, ... In the simplest case, if n = 2, then r = 6.Here we will use an equation of kind of Bernoulli as simplest equation.For the purpose we assume that c 0 = c 1 = c 3 = c 5 = 0; c 2 = p 2 ; c 4 = 2pq; c 6 = q 2 0, and search for a solution in the form: We note that at another choice of the parameters c j we can obtain another kind of dφ dξ , but we aim to use some simpler equation which has an exact solution according to the modified method of simplest equation.The substitution of (34) in (31) leads to the following system of relationships for the parameters of the solution (the system is of kind ( 4)): The system (35) implies that a 1 = 0 and q is a free parameter.The solution of (35) is: where The expression for the solitary wave depends on the solution of the differential equation in (34) and it is given by: for the case p > 0, q < 0 and for the case p < 0, q > 0, where a 0 , a 2 , p are given by (36).The values (37) or (38) describe 'kink' waves, which can be interpreted as fronts of some density propagation of the STAT5a signaling protein in the intracellular space, as we will show in the next paragraph of the paper.Two examples of these kinks are shown in Fig. 2 and Fig. 3.We note, however, that the figures give rather an illustrative idea of the possible spatial-temporal behavior of the considered object due to the lack of experimental data for the model parameters.

C. Discussions
In Section III we demonstrated that the lowest and highest steady state values of the STAT5a concentration are realized practically in the homogeneous case corresponding to the absence of appropriate response in the form of ERK or SOCS protein production initiated by the cell signaling.Then stabilization of the STAT5a signaling protein occurs.When, however, similar response (ERK or SOCS proteins increase) is available, the STAT5a concentration dramatically increases or decreases near the nucleus membrane and some inhomogeneous distribution of the STAT5a molecules takes place in the cytosol.This means that the diffusion involves in the process and it can be presented by (25) (or ( 26)).In this case the STAT5a concentration will propagate in the intracellular space as a "kink" wave, as we demonstrated in the previous paragraph.We choose the wave propagates along the positive direction of axis , i.e. -from the cell membrane to the nucleus one according to the direction of the signal transduction.Equation ( 25) (in particular (26) or (27) in view of the introduced traveling-wave ansatz) has the same steady state values for s as in the homogeneous case which we considered in Section III.However in inhomogeneous case, presented by (25) or (26), the lowest and the highest steady states are unstable, and the intermediate one is stable.Therefore the steady state values in inhomogeneous case change the type of their character in comparison with the homogeneous one.In this way at additional cell signals, initiating increase of the active ERK molecules (but in absence of additional SOCS protein production), the STAT5a concentration will increase from the lowest destabilized value s 0 1 to stabilized one s 0 3 , and density (concentration) jump will be observed.A similar example is shown in Fig. 2. On the other hand at additional new signal (ligand), which can inhibit the STAT phosphorylation (For example, the SOCS protein production increases, but there is not ERK protein production), the STAT5a concentration can decrease from its highest destabilized value s 0 2 to stabilized one s 0 3 , and density (concentration) drop will be realized in the intracellular space (Fig. 3).Thereby the inhomogeneous effects observed in the STAT5a protein dynamics are predetermined by the levels of the ERK and SOCS proteins.

V. Conclusion
In this paper, we have re-considered the reaction model of the cross talk between ERK and STAT signaling pathways.By applying a well-known QSSA theorem to the described model we have demonstrated that near the quasi-stationary state of the considered process its dynamics determines by the behavior of the phosphorylated STAT5a protein.Moreover, if the STAT5a concentration is homogeneous distributed in the cell cytoplasm its possible minimum or maximum levels will support stabilization of the cell signaling.But, if some inhomogeneous distribution of the STAT5a protein appears in the intracellular space as a result of external signals destabilization of the above mentioned levels will occur.In order to ensure the stabilization of the cell signaling process, the STAT5a density (concentration) will move in the cell from its unstable levels to stable ones.By applying the modified method of simplest equation to the corresponding reaction-diffusion model we have derived an analytic expression (solution) for the spatial STAT5a wave propagation .The 'kink' form of the considered protein density wave supposes eventual drop and jump concentration effects.These effects are a sign for cell inhomogeneity, which deviates the cell from the homeostatic norms and enables the initiation of pathological processes.