Reaction-Diffusion Modeling ERK- and STAT-Interaction Dynamics

The modeling of the dynamics of interaction between ERK and STAT signaling pathways in the cell needs to establish the biochemical diagram of the corresponding proteins interactions as well as the corresponding reaction-diffusion scheme. Starting from the verbal description available in the literature of the cross talk between the two pathways, a simple diagram of interaction between ERK and STAT5a proteins is chosen to write corresponding kinetic equations. The dynamics of interaction is modeled in a form of two-dimensional nonlinear dynamical system for ERK—and STAT5a —protein concentrations. Then the spatial modeling of the interaction is accomplished by introducing an appropriate diffusion-reaction scheme. The obtained system of partial differential equations is analyzed and it is argued that the possibility of Turing bifurcation is presented by loss of stability of the homogeneous steady state and forms dissipative structures in the ERK and STAT interaction process. In these terms, a possible scaffolding effect in the protein interaction is related to the process of stabilization and destabilization of the dissipative structures (pattern formation) inherent to the model of ERK and STAT cross talk.


INTRODUCTION
One of the features distinguishing a modern dynamics is its interest in framing important descriptions of the real processes in the form of dynamical systems. We call dynamical a system of first-order autonomous ordinary differential equations solved with respect to their derivatives. In some cases partial derivatives are included too and the corresponding systems are called spatial dynamical systems. The process of translation of observed data into a mathematical model in this case is called dynamical modeling (Beltrami [1]) and spatial one in particular. Dynamical systems belong to one of the main mathematical concepts. It is clear that dynamical systems constitute a particular case of the numerous mathematical models that can be built as a result of studies of the world that surrounds us. In view of the fact that there are different types of dynamical models, we restrict our considerations on none but models described by dynamical systems defined above.
The system analysis of intracellular processes and especially signaling, excitation, and mitosis (growth and division) in eukaryotes is so complex that it defies understanding by verbal arguments only. The insight into details of biochemical kinetics of cell functions requires mathematical modeling of the type practiced in the classical dynamics, that is, by dynamical systems. They are systems of differential equations arrived at in the process of studying a real phenomenon. In this paper we propose a dynamical modeling of intracellular processes. For this purpose the molecular mechanism of ERK (extracellular-signal-regulated kinase) and STAT (signal transducer and activator of transcription) pathways interaction is presented verbally and by a corresponding biochemical diagram. On this basis we write out a system of nonlinear ordinary differential equations (ODEs) expressing the kinetic mass action. Then we show at equilibrium that the ODEs become quadratic equations, whose solution describes the equilibrium concentrations. To understand how stable the equilibrium is, we use a small perturbation term to see how the differential equations governing the rate of change of the perturbation can be approximated. Next we use the standard Routh-Hurwitz condition to characterize the stability type of the steady state (equilibrium). What is of essential interest further is the question "how could we handle diffusion-reaction (partial differential) equations by first 2 EURASIP Journal on Bioinformatics and Systems Biology analyzing diffusion along one dimension, then proceeding to Turing bifurcation analysis?" We perform stability analysis on this reaction-diffusion system again by solving for the equilibrium and then studying its perturbations. At the end by analogy with the dynamical behavior of ERK and STAT interaction we propose a hypothetical scaffolding mechanism of the process.
The motivation and purposes above-mentioned lie in the following circumstance: on one hand the complexity of intracellular space is inscribed by the huge amount of interacting proteins and their molecular pathways and networks. On the other one, the heterogeneous distributions of protein concentrations in the form of cellular compartments play a crucial role in the regulation of all processes in the cell. In this way, cellular complexity is inherently space-temporal, described physically as reaction-diffusion processes not only between organelles and cytosol, but as a set of interactions between compartments and cytosol. The traditional approximation scheme of well-stirred reactor is a simplification due to the added complexity of modeling diffusion as well as the lack of straightforward experimental techniques to provide the necessary measurements needed to fully describe a spacetemporal model (Eungdamrong and Iyengar [2]). If the time resolution of the system is large enough, this approximation is valid for many materials with fast diffusion rates and/or small volumes. At this condition, diffusion acts simply as a mechanism to slow down the apparent associative or disassociative rate constant, and transport between compartments may be effectively treated as gradients between spatially averaged concentrations of the transported species. However, the concentration gradients of enzymes within cells that modulate signal transduction belie this simplification (Khurana et al. [3]; Holdaway-Clarke et al. [4]; Lam et al. [5]; Belenkaya et al. [6]). With experimental and technological advancements allowing finer temporal and spatial resolution, the development of space-temporal (i.e., reaction-diffusion) modeling intracellular kinetics to traditional systems biology has become much more tractable. That is why here we introduce both methodological foundation by proposing a specific technique of reaction-diffusion modeling and its computational implication to concrete example of ERK and STAT protein interaction. The specificity of this approach is also in the combining of an appropriate scheme of modeling with its analysis by the method of stability and bifurcation theory of dynamical systems. Similar approaches suggested that analyzing chemical systems were previously proposed in molecular chemistry (Lengyel and Epstein [7]). They obtained two-dimensional system of Turing type for the case of chlorine dioxide/iodine/malonic acid reaction and suggested hypothesis that a similar phenomenon may occur in some biological pattern formation process as it is in our case. In this sense our work could be considered as a confirmation of Lengyel and Epstein hypothesis. In a more general plan (ndimensional case) the problem of pattern formation is considered using rigorous mathematical terms in the paper of Alber et al. [8].
The approach in this paper takes into account the specificity of cell signaling of ERK-and STAT-pathways involved in a corresponding kinetic scheme different from those in the papers of Lengyel and Epstein [7] and Alber et al. [8] and applies appropriate mathematical methods (Lyapunov's stability and Tihonov's theorem). The significance and utility of our specific approach to modeling dynamically a possible scaffolding mechanism and dynamical nature of ERK and STAT interaction is discussed in the last two sections.

THE INTERACTION BETWEEN ERK AND STAT PATHWAYS: A DYNAMICAL MODEL
It is known that growth factors typically activate several signaling pathways. On this basis the specificity of biological responses is often achieved in a combinatorial fashion through the concerted interaction of signaling pathways (Pawson et al. [9]). The explanation is that many of the signaling pathways and regulatory systems in eukaryotic cells are controlled by proteins with multiple interaction domains that mediate specific protein-protein and protein-phospholipid interactions, and thereby determine the biological output of receptors for external and intrinsic signals. In the mentioned paper of Pawson et al. [9] the authors discuss the basic features of interaction domains, and suggest that rather simple binary interactions can be used in sophisticated ways to generate complex cellular responses. In the paper of Shuai [10], the protein STATs (signal transducer and activator of transcription) is 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 block access of STAT to receptors or inhibit JAKs or both (Alexander [11]). (JAK is an abbreviation of Janus kinase). On the other hand, SOCS-3 can bind to and sequester such named Ras-GAP (Cacalano et al. [12]). The suppressors of cytokine signaling (SOCS, also known as CIS and SSI) are encoded. By immediate early genes they act in a feedback loop to inhibit cytokine responses and activation of signal transducer and activator of transcription (STAT). The activity of signal transducer activator of transcription 5 (STAT5) is induced by an overabundance of cytokines and growth factors and resulting in a transcriptional activation of target genes (Buitenhuis et al. [13]). STAT5 plays an important role in a variety of cellular processes as immune response, proliferation, differentiation, apoptosis. What is of interest from medical point of view, aberrant regulation of STAT5 has been observed in patients with solid tumors, chronic and acute myeloid leukemia.
In the papers of Wood et al. [14]; Pircher et al. [15], it is suggested that the STAT5 functional capacity of binding DNA could be influenced by the mitogen-activated protein kinase (MAPK)-pathway. Moreover, it is known that the serine phosphorylation of signal transducers and activators of transcription (STAT) 1 and 3 modulates their DNA-binding capacity and transcriptional activity. In a later Nikola Georgiev et al. 3 paper of Pircher et al. [15] the interactions between STAT5a and the MAPKs (extracellular signal-regulated kinases ERK1 and 2) are analyzed. In vitro phosphorilation of the glutathione-S-transferase-fusion proteins using active ERK only worked when the fusion protein contained wild-type STAT5a sequence. Transfection experiments with COS cells showed that kinase-inactive ERK1 decreased GH stimulation of STAT5-regulated reporter gene expression. These observations show for the first time a direct physical interaction between ERK and STAT pathways. They identify also serine 780 as a target for ERK.
From the results described in the work of Pircher et al. [15] a model for interaction between ERK and STAT5a in CHOA cells can be derived (Figure 1), we call it a model of Pircher-Petersen-Gustafson-Haldosen or PPGH-model (diagram). As it is seen from Figure 1, in unstimulated cells STAT5a is complexed with inactive ERK that binds to STAT5a via its C-terminal substrate recognition domain to an unknown region on STAT5a. Then via its active site it binds to the C-terminal ERK recognition sequence in STAT5a. On the other hand, upon GH stimulation, MEK activates ERK through phosphorilation of specific threonine and tyrosine residues in ERK. As shown in the paper of Pircher et al. [15], the cytosol and nuclear extracts of in vitro cells were analyzed using Western blotting method; by using antibodies against ERK1/2, active ERK1/2, and STAT5a. The relation in Figure 1 was derived from the Western blotting qualitative results. Later, other publication revealed the insides of the two ERK/MAPK and JAK/STAT pathways. It is already known that during growth factor stimulation, the ERK phosphorylation cascade is linked to cell surface receptor tyrosine kinases (RTKs) and other upstream signaling proteins with oncogenic potential (Blume-Jensen and Hunter [16]). The MAP kinases ERK1 and ERK2 are 44-and 42-kDa Ser/Thr kinases, with ERK2 levels higher than ERK1 (Boulton et al. [17,18]).
From the diagram in Figure 1 we can write the following system of ordinary differential equations for the kinetics of STAT5a/S phosphorylation and ERK activation, described by concentration variables e 1 , e 2 , s 1 , s 2 denoting concentrations of ERK-inactive, ERK-active, STAT-and STATphosphorylated, respectively. It has the form where k 1 is proportional to the frequency of collisions of ERK and STAT protein molecules and present rate constant of reactions of associations; k 2 and k 3 are constants of exponential growths and disintegrations; I > 0 inhibitor source respectively. The source I inhibits the phosphorylation of non- proteins. Thus mathematically, as a first approximation we can write where Σ is a constant concentration of SOCS proteins and k is a reaction rate constant of inhibition, respectively. It is clear that if Σ increases, the term I increases too and vice versa.
To analyze (1) we pay firstly attention that only two equations of the four ones are independent. It is easy to show that between the concentrations e 1 , e 2 , s 1 , s 2 there exist the relations where are initial values in the interval (0,1) of the sums of corresponding concentrations of inactive and active ERKs and nonphosphorylated and phosphorylated STATs. The relations present the steady state of (1). The notation e in (4)-(5) is a noninteracting part of the concentration of ERK proteins. Moreover, s 0 2 is a positive real root of the quadratic equation where

EURASIP Journal on Bioinformatics and Systems Biology
The eventual negative or complex roots have not physical sense. From the expressions (7) for β and γ we conclude that they become respectively positive and negative with large absolute values when Σ is large. Then, from the formula of the roots of (6) s 0 it follows that in this case (Σ is sufficiently large) (s 0 2 ) 1 is always positive and (s 0 2 ) 2 is negative. Moreover we can choose (s 0 2 ) 1 large by choosing corresponding large Σ (high concentration of SOCS proteins). We could do all this independently of the values of E and S (e.g., E sufficiently small and S large). The smallness of E follows from the consideration that the inactive ERK concentration could in principle contain both participating e 1 and not participating e parts in the ERK and STAT interaction.
Further we replace e 1 and s 1 from (3), respectively, in the second and fourth equations of (1). As a result we obtain the two-dimensional system having a steady state It is clear that if the equilibrium (10) of the twodimensional system (9) is stable, then the equilibrium (5) of the four-dimensional system (1) is stable too. In order to analyze the stability of the equilibrium (10) we linearize (9) around (10) by substituting the changes where ξ, η are variations (disturbances) around the steady state. Then the variation equations of the model (9) take the form where for the coefficients in the right-hand side, the following formulas are valid The Routh-Hurwitz conditions for stability of the steady state (10) have the form In view of the first formula of (10) we can conclude the following.
(1) At the absence of noninteracting ERK proteins, when E = e 1 + e 2 is strictly valid, the conditions (14) are satisfied, because in this case the inequalities E e 0 2 > 0, S s 0 2 > 0 always hold and the coefficients k 1 , k 2 , k 3 are positive too (by definition).
(2) When the concentration of noninteracting ERK proteins is sufficiently large, the inequalities (14) become opposite.
(3) For small E, large S and Σ, the following relations are possible: a > 0, c > 0, b < 0, d < 0 under condition that (14) hold. These are necessary conditions for such named Turing bifurcation of the distributed version of the model (12).
If the disturbances ξ and η are sufficiently small, then the system (12) can be reduced to the following linear oscillator with attenuation and under external influence where the new variable x(t) presents both signals ξ and η.
The function f (t) presents some permanent external influence on ξ and η. The analysis of (15) is well known and we present here only the most essential of the results. The functions f (t) and x(t) can be presented in the form of the following Fourier-integrals: where the functions F(ω) and X(ω) are spectral densities of the functions f (t) and x(t), respectively. By substituting (16) in (15) we obtain from where we find If the attenuation γ is small, what seems possible in view of the formulas (14), then X(ω) can be too large, when the external frequency ω is near the resonant frequency ω 0 . Thus in the Fourier spectral density of x(t) the most large are X(ω 0 ) and X( ω 0 ), when we can talk about resonance phenomenon in signaling.

MULTICOMPONENT ONE-DIMENSIONAL DYNAMICAL SYSTEM WITH DISTRIBUTED VARIABLES
The role of diffusion in reaction-diffusion systems of the cell becomes significant when reactions are relatively faster Nikola Georgiev et al.

5
(but not too very) than diffusion rates and is known in the literature as spatial distributed process. Sometimes the term crowding is used to denote a more specific type of spatial distribution (Takahashi et al. [19]). The physicochemical essence of this phenomenon lies in the circumstance that the state of phosphorylation of target molecules with spatially separated membrane-localized protein kinases and cytosolic phosphatases depends essentially on diffusion (Kholodenko et al. [20]). The crucial coupling of diffusion and noise is implied by the fact that subcompartments diffusively formed by localized proteins can definitely alter the effect of noise on signaling outcomes (Bhalla [21]). The very high protein density in the intracellular space, commonly called molecular crowding, can augment the spatial effect. Consequently, molecular crowding can also alter protein activities and break down classical reaction kinetics (Schnell and Turner [22]). In the remainder of this article, we develop a mathematical approach that can be used to model and simulate the consequences of spatial distribution. Although we will only consider MEK/ERK and JAK/STAT-signaling pathways, most discussions in this paper should also be applicable to other intracellular phenomena. They involve reactiondiffusion processes as EGF signaling pathway, interleukins IL2, IL3, and IL6 signaling pathways, inhibition of cellular proliferation in Gleevec, PDGF signaling pathway, or TPO signaling pathway. It is known that signalling pathway MEK/ERK can be activated and regulated by dynamic changes in their organization both in time and space. The JAK-STAT signaling cascade is also characterized by the activation of a JAK-kinase that is bound to the cytoplasmic domain of a cell surface receptor such as the erythropoietin receptor (EpoR) (Swameye et al. [23]). Moreover, in the paper of Ketteler et al. [24] it is shown that a receptor harbouring the GFP (Green Fluorescent Protein) inserted near the two STAT5 binding sites in the EpoR cytoplasmic domain retains full biological activity. In a similar way, we know from Kolch [25] that the ERK pathway features dynamic subcellular redistributions closely related to its function. As a rule the activation of Raf-1 and B-raf ensue with the binding to Ras resulting in the translocation of Raf from the cytosol to the cell membrane. Many questions arise however in both JAK/STAT and MEK/ERK for clarifying dynamic details of time-space effects. In order to answer them we should develop a general approach to modeling the special relocalization process in the cell.
The variation of signal components along time and space (in the cell) can be described by such a named diffusionreaction equation, having the form where c is the concentration of the signal component (as a rule-some protein), t is the time, k is a diffusion coefficient of signal molecules, f is a velocity of production and consumption of the signal component, what is in principle nonlinear function of c (Georgiev et al. [26]). In this way (*) is a nonlinear differential equation in partial derivatives. Its deduction can be found in the book of Berg [27].
The diffusive coefficient predetermines the range of diffusion signal components by the well-known formula for the dependence of the range radius on the squared root of the diffusive coefficient. It is known that the signal network participating in the morphogenesis of the biological development is considered as dependent on the local activation of the components and their global inhibition (Berg [27]; Nagorcka and Mooney [37]; Painter et al. [28]). What is of interest in our paper is the possibility that similar space localized reactions can be modeled by small diffusive coefficients for the components with positive feedback loops (activation) and by large diffusion coefficients for the components with negative feedback loops (inhibition). Concerning these, here the concepts of stability and instability are widely treated in general sense and applied to corresponding ERK and STAT spatial models. For this purpose, Lyapunov's method of first approximation is systematically applied. In the literature, the stability analysis of reaction-diffusion equations (rde) is often connected with the realization of possibility that dynamical systems in the infinite phase space are to be reduced to lowdimensional systems. These are problems of reduction possibly solvable by such named methods of projection, based on the known Fredholm theorem (Iooss and Joseph [38]).
In this section we introduce a generalization of the monocomponent rde in the form (*) to multicomponent case of many concentrations. For this purpose we define firstly some general notions. We call systems with distributed variables when the connections between neighbor points of space are taken into account by the diffusion law of molecular motion from the higher to lower concentrations. In onedimensional case (not monocomponent) when the diffusion occurs along space coordinates, the full system of differential equations by accounting the diffusive terms can be written in the form where the functions Q i (x) define dependence of the concentrations c 1 , c 2 , . . . , c n on the space coordinate x, and the nonlinear functions f i (c 1 , c 2 , . . . , c n ) in the right-hand side correspond to a "point" model, that is, with concentrated parameters. The very spatial distribution in the cell is presented by reaction-diffusion process of interaction between proteins and protein-complexes of the signaling pathway and takes place in some intracellular volume described below. Let us assume that the solution of (19) has the form In order to find in explicit form the functions Q i (x), we consider the signal pathway as being contained in a simple intracellular domain having the form of long narrow tube with a length l and cross section S (Figure 2). In this tube we separate an elementary volume ΔV with limit coordinates x and x + Δx.
where D is a diffusion coefficient, defined by the properties of solution substances. In spite of the other limit of the volume with coordinate x + Δx, in the opposite direction and during the same time interval, it diffuses a mass In this way, the total mass variation in the elementary volume ΔV at the expend of diffusion is and the variation of concentration c i is presented by By limit transition to Δx 0 we obtain By definition, in the absence of biochemical reactions in correspondence with (19) we have Q i = lim(Δc i /Δt), when the limit transition Δt 0 takes place. Thus, at the same transition we can write where the quantities Q i have the same physical sense as in (1). Therefore, the distributed system (1) in case of onedimensional diffusion has the form where the nonlinear functions f i (c 1 , c 2 , . . . , c n ) correspond as before to the point model and D i (∂ 2 c i (x, t)/∂x 2 ) correspond to the diffusion transport between the neighbor volumes. Equation (27) presents a system of nonlinear differential equations in partial derivatives. In order to analyze qualitatively and solve quantitatively these equations, it is necessary to fix some initial conditions in the form of initial distribution of the unknown concentrations c i along the space coordinate x in the moment t = 0, that is, Moreover, the values of concentrations at the boundary of the reaction volume V of the signal pathway must be given too. If the reaction volume of the pathway is sufficiently large, then it is not necessary to take boundary conditions. It is of interest to know the cases when (27) can be reduced to a system with concentrated parameters (point system). They are the following.
(1) When all coefficients of diffusion vanish, that is, D i = 0. In this case the protein molecules and protein complexes will not collide each the other and the biochemical reactions of the signaling pathway will not occur. A signal pathway does not exist.
(2) If the diffusion coefficients are very large (D i ½), the diffusion velocity will be large with respect to the rate of biochemical reactions. Then before the essential variation of concentrations at the expense of the biochemical reactions, the protein molecules and protein complexes will displace through the whole pathway volume. Thus after some very short time of relaxation, the solution of (27) will approach very near to the solution of corresponding model with distributed variable of the pathway.
(3) When the outer conditions (out of the reaction volume of the pathway) and initial conditions are homogeneous in whole volume, that is, ϕ i (x) = ϕ i = const, that means the diffusion is absent and it is also sufficient to consider only a point system (with concentrated variables).
The specific applications of systems with distributed variables (concentrations) of type (27) to mathematical descriptions of spatial relocalization processes in cell present difficult problems. That is why we will consider only some simple examples in order to illustrate the application of similar systems to describing intracellular processes. For this purpose we should note first of all that the biological systems with distributed variables (including also signaling pathways) belongs to such called active distributed systems. They are characterized by a sequence of properties called qualitative particularities. These are the emergence of nerve excitation (action potential) in the nerve cell, autocontraction of the cardiac cell and other instabilities and bifurcations, leading to various regimes of functioning in cell differentiation and proliferation. It is reasonable to expect that paradigmatic models of type of (27) can be used to describing processes of proteins distribution in the cell at signaling pathway level. What is of interest in this case is that the form of nonlinear functions f i , the relationships between parameters and their values determine the regime of system functioning: stable, not depending on time, nonhomogenous space solutions, traveling impulses, synchronic self-oscillations of the whole pathway or of the separated parts only.
In the next sections we will restrict the consideration only to the following basic stages of analyzing distributed systems (27).
(1) Finding steady state homogeneous or nonhomogeneous in the space solutions constant along the time.
(2) Studying the stability of the found steady state solutions.
(3) Evolution of distributed system along time and appearance of dissipative structures in the signaling pathway.

STABILITY ANALYSIS OF THE HOMOGENEOUS STEADY STATE OF ERK AND STAT INTERACTION WITH DIFFUSION
Now we apply the procedure developed in the previous section to the dynamical model of ERK and STAT interaction in the form (12). As a result we obtain the following twodimensional system with distributed parameters: where r is the space coordinate from the cell membrane to the nucleus, D ξ , D η are coefficients of diffusion of the concentration deviations (disturbances) ξ, η respectively. In order to analyze qualitatively and solve quantitatively these equations, it is necessary to fix some boundary conditions for the gradients of concentrations at the cell membrane and nucleus in the form ∂ξ ∂r r=0 r=l = 0, ∂η ∂r r=0 where l is the distance between the membrane and nucleus. The steady state of (29) is homogeneous and has the form It is equivalent to the homogeneous equilibrium e 0 2 (t, r) = k 3 s 0 2 + kΣ /k 2 , s 0 2 (t, r) = s 0 2 .
of the model reaction-diffusion system of ERK and STAT interaction Here D s = D ξ and D e = D η . Our mathematical model (1), (2), (29)-(33) of reactiondiffusion is based on the oversimplified model of (Pircher et al. [15]), Figure 1, obtained by qualitative and not quantitative Western blotting. That means the mathematical analysis of our model despite of its high complexity (e.g., high number of parameters of the system) must be also qualitative and not quantitative one. In view of this we will use the language of the nonlinear dynamical systems theory, which is qualitative and very similar to the traditional biochemical one, being verbal and needing mathematical accuracy (in qualitative sense). Certainly, similar approach requires verification of a qualitative correspondence between the effects of theoretical predictions and experimental measurements. In particular, it is in concordance with claiming qualitative scaling relationship in terms of Tichonov's theorem (Tichonov [29]), as we will do further.
To investigate the stability of (31), (32), we should obtain the solutions of the linear system which is valid for small disturbances ξ, η. If the solution of (34) attenuates, then the homogeneous steady state (31) (or (32)) is stable. Otherwise, it is unstable and an emergence of dissipative structures is possible in principle.
For infinite one-dimensional space the value of the wavelength λ changes continuously from 0 to ½, and in case of segment (as it is our case), λ takes discrete values. The complex frequency p is defined from the quadratic equation Consider a relationship between the real part of roots of (36) and the parameter u = (2π/λ) 2 (the square of wave number). Let us now accept that D ξ > D η , where D ξ is the diffusion coefficient of the molecules of STAT5a protein, which are larger than that of ERK noted by D η . If bc > 0, then both roots p 1,2 are real numbers for every value of λ (Figure 3). If bc < 0, then p 1,2 are complexly conjugated numbers for wavelength in the interval 4π 2 /u 4 λ 2 4π 2 /u 3 (Figure 3), where u 3 and u 4 are equal to In every graph of Figure 3, we can separate 3 regions: (I) both roots p 1,2 have positive real part, that is, Re p 1,2 > 0; (II) one of the roots has a positive and the other-negative real parts, that is, Re p 1 > 0, Re p 2 < 0; (III) both roots p 1,2 have negative real parts, that is, Re p 1,2 < 0. By using the terminology of the qualitative theory of dynamical systems, we say that the linear system (34)  on the straight line of the parameter u are defined by the values of u 1 , u 2 for which one of the real parts Re p 1,2 becomes zero: It can be shown that under perturbations with wavelength of the region I in nonlinear distributed system can emerge waves with final amplitude; under perturbations with wavelength of the region (II) spatially periodical steady states regimes (such named dissipative structures) emerge.

STABILITY ANALYSIS OF THE INHOMOGENEOUS STEADY STATE OF ERK AND STAT INTERACTION WITH DIFFUSION
Let us consider again the distributed nonlinear model (29) of ERK and STAT interaction under boundary conditions (30) from the previous section. We pay attention that ξ and η are finite deviations (disturbances) of the STAT and ERK protein concentrations from the steady state values (10). The last obtained by equating to zero the right-hand sides of the ERK and STAT interaction model with concentrated parameters (i.e., ordinary differential equations) (9). There k 1 is considered as a relatively small (with respect to a ) coefficient, proportional to the frequency of collisions of ERK and STAT protein molecules and presents a rate constant of reactions of associations, and Σ is sufficiently high (to assure s 0 2 > 0). In steady state approximation the model (29) takes the form Further we assume the inequality D ξ D η in view of the fact that the ERK molecule is smaller than STAT one (Pircher et al. [15,31]). This circumstance can be related to the fact that STAT pathway tends to be much more rapid than the ERK one. For this purpose we consider the first equation of (39) to be linear with respect to ξ and can be treated as an attached system in accordance with the Tichonov's theorem (Tichonov [29]). Next we take into consideration that η is a sufficiently small "constant." Thus the attached system has a stable steady state of the center type (then well-known Lyapunov's definition of stability is satisfied). After replacing the steady state value of ξ from the first equation in the second one (the degenerate system), the last is obtained in the form The corresponding reaction-diffusion equation is under the same boundary conditions (42). This cubic polynomial approximation means we accept a weak nonlinearity (but not linearization) of the model (29), that is, k 1 is sufficiently smaller than a or k 2 , k 3 to assure the approximation validity. The last inequalities follow from the biochemical consideration that the processes of ERK inactivation and STAT dephosphorylation are faster than that of ERK and STAT interaction. The last being of molecular recognition type (Pircher et al. [15,31]) with possible scaffolding mechanism to be assumed further.
Let us now substitute in (43) the perturbation solution η(t, r) = η(r) + ω(t, r), where η(r) is an inhomogeneous steady state solution of (41) and ω(t, r) is a small variation (perturbation). We obtain the next variation equation under initial condition (playing the role of a dissipative structure in this case) and boundary conditions ∂ω ∂r r=0 By applying the standard procedure similar to that in the previous Section 4, the solution of (44) can be obtained in the form ω(t, r) = ½ n=0 a n e Q(η)t cos λ n r, where a n = 2 l l 0 ϕ(r) cos λ n r dr, λ n = nπ l 2 , Next we denote where θ, τ, γ are positive numbers in view of the relations being valid at the absence of noninteracting ERK proteins. Then the expression (49) takes the form Here are the roots of the quadratic polynomial Q(η). There are two negative steady state values of the deviation η from the steady state value of the concentration e 0 2 assumed to be larger than the corresponding deviations. It is easy to show that Q(η) is negative when the steady state concentration is out of the interval between the two roots mentioned. In this case the perturbation solution (47) attenuates and the dissipative structure (45) is stable, thus it could really exist. For a steady state deviation smaller than the bigger root and larger than the smaller one, Q(η) is positive if the structural wave number λ n or diffusion coefficient D η is sufficiently small and then the dissipative structure (45) is unstable and disappears. Thus too low and too high steady state concentrations are indicative for the dissipative structures existence, but the average ones are not. Following this, in the next section it will be shown how the well-known scaffolding effect (Stewart et al. [32]; Schaeffer et al. [33]; Teis et al. [34]) can be related to the described behavior of η (activated ERK).

HYPOTHETICAL MECHANISM OF STAT SCAFFOLDING ERK PATHWAY
In terms of the above described stability analysis, the magnitude of initial disturbance η of activated ERK depends critically on its own value. Corresponding initial values of η can amplify or attenuate in a regime of instability or stability respectively. The dynamical interpretation of ERK criticality consists in the effect above theoretically established that in some interval of ERK concentration the ERK pathway is unstable. That means initial concentration of ERK belonging to this interval does not conserve its amplitude but amplifies. Thus ERK pathway is sensitive at intermediate concentrations of ERK. Out of this interval of average ERK concentrations, the ERK pathway becomes insensitive, that is, the distribution conserves its magnitude unchanged. This purely qualitative consequence from our model can be explained physically by hypothetical STAT scaffolding mechanism of ERK signaling, presented in Figures 4 and 5. Before the latter mechanism can be addressed, we need to define a scaffold as a protein whose main function is to bring other protein together for them to interact. Such a protein usually has many protein binding domains what is not yet established for STAT. In the basic work (Pircher et al. [15]) it is mentioned about "unknown region on STAT5a." Concerning that we accept the hypothesis that STAT may has several  binding domains and we try to draw some conclusions from this assumption. Papers by Bray and Lay [35] and Levchenko et al. [36] have provided corresponding insights in general sense into this hypothesis through computer simulations of signaling pathways with scaffolds. On the basis of these studies the first idea we can relate to STAT scaffolding mechanism is presented in Figure 4. It illustrates the principle of balance: adding too much ERK concentration we can decrease the output of ERK scaffolded cascade, just as adding too much scaffold STAT can ( Figure 5). The analogy of the presented simple mechanism with the dynamical behavior of ERK signaling is evident: in both cases ERK pathway amplifies signal for intermediate concentration of scaffold STAT and does not amplify it for low and high concentrations. In Figure 5 it is seen again a scheme like combinatorial inhibition. Signaling down scaffolded ERK cascade is a question of balance: if there is too small STAT concentration, ERK signaling will be low (left). At an intermediate STAT concentration, the ERK signaling will be high (center). Once the STAT concentration exceeds that of the ERK it binds, the signaling begins to decrease (right).
The most important question now is whether ERK and STAT interaction really exhibits the scaffold mechanism predicted in this section. With the exception of unknown number of binding domains of scaffold STAT, for which necessity to be measured there is good experimental evidence, there is not much principal objection against the hypothetical mechanism suggested here.

CONCLUSION
The present analysis shows that diffusion (together with corresponding biochemical reactions) is likely to play a critical role in governing the space-temporal behavior of ERK and STAT interaction system and should not be ignored. In terms of the reaction-diffusion interaction in ERK and STAT dynamical model presented here, the effect of protein scaffolding can be related to a destabilization of inhomogeneous distributions of protein concentrations.
In view of the fact that the modeling parameters are usually gathered from biochemical experiments on purified components while functional effects arise from cell physiological experiments, one does not aim at numerical agreement between experimental data of scaffolding effect and some modeling prediction. Instead, the modeler should aim for correct "scaling relationship" in qualitative sense (relatively large and small). The gradual refinement of a corresponding dynamical model (e.g., (29)) should be an iterative process. Once the time scale of activation and deactivation processes of ERK and STAT interaction are qualitatively established, one can investigate the robustness (including possible stabilities and instabilities of steady states, oscillations and spatial patterns of the model) by bifurcation analysis of system behavior in dependence on its parameters. The use of qualitative analysis (in terms of qualitative theory of dynamical systems) would be in a possible confirmation or rejection the structural stability of the signaling pathway under consideration. It seems established that the most signaling networks of protein interactions in the cell are structurally stable within an order of magnitude variation in kinetics parameters. That is, our model should not "blow up," (i.e., concentration of a reactant goes to infinity or zero) when the parameters are increased or decreased by less than ten-fold. If, on the other hand, the structural stability analysis predicts that some components (concentrations) are sensitive to small perturbations in model parameters (rate constants), this would not necessarily mean the model is incorrect. In similar case it could suggest that the components might be susceptible to external perturbation, and experimental means to manipulate the biochemical activities of the molecule should be used to verify theoretical prediction.
The present analysis shows that diffusion is an essential part of cellular complexity inherent to space-temporal description not only as a set of complex protein networks within organelles and the cytosol, but also as that of interactions between compartments and the cytosol. A new generation of microscopic techniques capable of resolving the intracellular localization of proteins provides evidence of the biological significance of process dynamically evolving in both space and time. Computer simulation of physicsbased models, coupled with quantitative space-temporal data will allow cell biologists to rigorously develop and test complex hypothesis of the dynamical nature of signaling pathways.

ACKNOWLEDGMENT
This work was supported by the European Community as part of the FP6 COSBICS Project (512060).