- Research article
- Open Access
- Published:

# A systems biology approach to analyse amplification in the JAK2-STAT5 signalling pathway

*BMC Systems Biology*
**volume 2**, Article number: 38 (2008)

## Abstract

### Background

The amplification of signals, defined as an increase in the intensity of a signal through networks of intracellular reactions, is considered one of the essential properties in many cell signalling pathways. Despite of the apparent importance of signal amplification, there have been few attempts to formalise this concept.

### Results

In this work we investigate the amplification and responsiveness of the JAK2-STAT5 pathway using a kinetic model. The recruitment of EpoR to the plasma membrane, activation by Epo, and deactivation of the EpoR/JAK2 complex are considered as well as the activation and nucleocytoplasmic shuttling of STAT5. Using qualitative biological knowledge, we first establish the structure of a general power-law model. We then generate a family of models from which we select suitable candidates. The parameter values of the model are estimated from experimental quantitative time-course data. The final model, whether it is conventional model with fixed predefined integer kinetic orders or a model with variable non-integer kinetic orders, is selected on the basis of a good agreement between simulations and the experimental data. The model is used to analyse the responsiveness and amplification properties of the pathway with sustained, transient, and oscillatory stimulation.

### Conclusion

The selected kinetic model predicts that the system acts as an amplifier with maximum amplification and sensitivity for input signals whose intensity match physiological values for Epo concentration and with duration in the range of one to 100 minutes. The response of the system reaches saturation for more intense and longer stimulation with Epo. We hypothesise that these properties of the system directly relate to the saturation of Epo receptor activation, its low recruitment to the plasma membrane and intense deactivation as predicted by the model.

## Background

Cellular signal transduction is accomplished by networks of interacting proteins that detect, modulate and transfer cellular signals which control gene expression. A prime example of this is tumour progression and certain oncogenic processes, which directly relate to dysfunctions in signal transduction networks [1, 2]. So far the use of mathematical modelling in cell signalling has been limited by the availability of suitable experimental data. However, the systematic development of experimental techniques enabling the generation of time-resolved quantitative data [3–5] facilitates the identification of dynamic pathway models and their parameter values by fitting them to experimental time course data.

The amplification of signals is considered one of the essential properties in most of the cell signalling pathways [6]. The notion of amplification as an increase in the intensity of the signal through the signalling cascade is generally accepted and used to characterise such systems. Surprisingly, there is little work in which a formal definition of amplification is proposed and used together with mathematical models to analyse signalling systems [7–10].

The Janus kinase – signal transducer and activator of transcription (JAK-STAT) pathways are one of the best-studied cell signalling pathways [11, 12]. The JAK2-STAT5 pathway is activated through various receptors, including the erythropoietin receptor (EpoR). Cytokine-activated phosphorylation of EpoR is mediated by the cytosolic kinase JAK2 which is associated with the cytoplasmic domain of EpoR. Upon binding of the hormone erythropoietin (Epo), JAK2 is activated and phosphorylates several tyrosine residues within the cytoplasmic domain of EpoR [13]. Subsequently, the transcription factor STAT5 is recruited to the activated receptor, becomes phosphorylated and thereby gets activated. Upon activation STAT5 homodimerises and migrates to the nucleus, where it initiates the transcription of target genes (see Figure 1 for a simplified representation of this pathway).

Previous work on data-based mathematical modelling of the core module of the JAK2-STAT5 signalling pathway [14] revealed nucleocytoplasmic cycling as an essential building principle of this pathway to closely couple gene transcription to receptor activation. Sensitivity analysis of the model showed that, surprisingly, not the first step of the pathway, i.e. phosphorylation of STAT5 at the receptor, but the parameters describing the shuttling through the nucleus have the major influence on transcription. The prediction of the outcome of an independent experiment based on this theoretical finding could be confirmed. Analysis of the model variables, i.e. unphosphorylated and phosphorylated monomeric STAT5 in the cytoplasm, and dimeric STAT5 in the cytoplasm and the nucleus, demonstrated that in a first round of activation nearly all accessible STAT5 molecules are phosphorylated. Thereby, it was shown that nucleocytoplasmic shuttling serves as a recycling step for the limited pool of STAT5 molecules, thus identifying an implementation of the strategies used by the cell to save energy and resources. The present paper extends this previous work by including a simplified description of the dynamics of EpoR at the plasma membrane (recruitment, degradation, and deactivation) and by analysing amplification and responsiveness of the pathway.

As an alternative to models that are derived on the basis of conventional mass action kinetics with predefined fixed integer kinetic orders, we use a more general power-law formulation that allows for non-integer kinetic orders [15] with the following structure:

Where *X*_{
i
}represents any of the *n*_{
d
}dependent variables of the model (e.g., proteins or phosphoprotein concentrations or levels of gene expression). Here, the biochemical rate *j* is expanded as a product of a rate constant (*γ*_{
j
}) and the *p* variables of the system to characteristic kinetic orders (*g*_{
jk
}), while *c*_{
ij
}are the stoichiometric coefficients of the system describing mass conservation. The main difference between power-law models and conventional ODEs models used in systems biology is that kinetic orders can have non-integer values. There are two main reasons to allow non-integer kinetic orders: firstly, reactions in non-homogenous environments lead to non-integer kinetic orders [16–18] and secondly in the absence of data on the detailed reaction mechanisms one is often forced to condense several steps into simplified representations. This aggregation of information is conveniently represented by power-law expressions, although other alternatives are possible [19, 20]. In power-law models, the kinetic orders are parameters of the model and must be estimated from experimental data. Negative values for the kinetic order represent inhibition, while a zero indicates that the variable does not affect the described process. When positive values are considered for a kinetic order, several alternatives are possible: values between zero and one represent a saturation-like behaviour for the rate modelled, and with values higher than one the rate equation models cooperative processes. A kinetic order equal to one means that the system behaves like conventional mass-action kinetic model. By allowing non-integer, positive or negative, kinetic orders we consider for the same model structure a larger class of kinetic models from which we can select a suitable candidate.

However, the increased generality comes at a price when it comes to the estimation of parameter values from experimental data for larger networks. While for conventional models the kinetic orders are decided *a priori*, in power-law representations the number of parameters that must be extracted from data is therefore larger. The identifiability problem means that there could be multiple sets of parameter values for which the model fits the data equally well. Analytical approaches for the inference of identifiability are limited to rather low dimensional systems [21], while approaches that are applicable to large systems [22, 23] rely on linearisation that question their reliability. Recent results based on non-parametric bootstrap-based approaches [24] will allow for reliable identifiability analysis also for high dimensional systems.

Next, we discuss the structure and equations of the mathematical model proposed and the process by which parameter values were estimated. In addition, we performed a model-based analysis of the responsiveness and amplification of the system, for which we propose a formal definition of these variables. The experimental and computational methods used to set up the model for the JAK2-STAT5 pathway are summarised at the end of the manuscript.

## Results and Discussion

### Mathematical Model and Parameter Estimation

In our model, EpoR and JAK2 were assumed to form a stable complex, EpoR/JAK2, for all biochemical processes included. All variables in the model describing the considered states of EpoR refer to populations of the receptor at the plasma membrane. Based on experimental data obtained during our investigation (data not shown), we assume in our model that the dynamics of the intermediate state of the receptor in which only JAK2 is activated is negligible in the description of the system. Thus, two possible states for the EpoR/JAK2 complex were considered: i) EpoR/JAK2 not bound to Epo, *EJ*; and ii) activated Epo-bound EpoR/JAK2 complex, *pEpJ*. The receptor activation was modelled with a term depending on the Epo concentration and the amount of non-activated receptor complex at the plasma membrane $\left({\gamma}_{2}\cdot E{J}^{{g}_{1}}\cdot Ep{o}^{{g}_{2}}\right)$ A rate equation accounting for the degradation of activated Epo receptor was also included in the model. Since the *k*_{
off
}of the Epo binding to the human EpoR is much smaller than the *k*_{
on
}[25], we assumed that effects of the release and reactivation of the related murine EpoR can be neglected, which was also tested in preliminary models. Actually, initial explorations considering the reactivation of EpoR were not consistent with the experimental data used in this work (data not shown). The structure of the model was completed by the inclusion of terms that describe the recruitment of new EpoR/JAK2 complexes to the plasma membrane and the degradation of non-activated EpoR/JAK2 complexes.

Three possible states have been considered for STAT5 in the model: i) non-activated and monomeric STAT5 in the cytosol, *S*; ii) activated STAT5 in the cytosol, *DpS*; and iii) activated STAT5 in the nucleus, *DpS*_{
nc
}. STAT5 activation comes through phosphorylation and dimerisation. Our preliminary results suggested that dimerisation is a fast processs and therefore the slower phosphorylation of STAT5 in the receptor complex leads the activation process. Thus, model reduction was applied and the aggregation of both processes was then represented with a unique rate equation (${\gamma}_{5}\cdot {S}^{{g}_{6}}\cdot pEp{J}^{{g}_{7}}$, see Additional file 1). Since the dimerisation of STAT5 is considered a fast process, the model assumes that the protein dimerises immediately after monomeric phosphorylation and the variable that describes monomeric phosphorylated STAT5 is neglected. We note that STAT5 can only be considered activated after the dimerisation of two phosphorylated monomers. Experimental data describing the translocation of activated STAT5 to the nucleus and its dynamical behaviour inside the nucleus are currently not available, and therefore the differential equations describing such processes have not been formulated in detail. Thus, we model the fraction of STAT5 inside the nucleus with a single state variable. The processes considered in the model are the activation of STAT5 by the activated receptor complex EpoR/JAK2, the translocation of cytosolic activated STAT5 to the nucleus, and the deactivation and subsequent translocation of nuclear STAT5 back to the cytoplasm. In accordance with [14], only nuclear activated STAT5 can be deactivated in our model. The total amount of STAT5 is supposed to be constant.

Finally, the concentration of Epo in the extracellular medium, *Epo*, is considered the input signal of the system. The resulting model has five dependent variables and one input variable. Figure 1 illustrates the structure of the model proposed. Only the states of the proteins and the processes that have been discussed are included. The differential equations for the model were formulated using power-law terms:

In this model we assume that, after the long period of starvation before starting the experiment, the system reaches a steady-state level of EpoR at the plasma membrane. Under these conditions and in the absence of stimulation, there is equilibrium between receptor recruitment and degradation, which makes the rate constants in both terms (recruitment rate and first order degradation for the receptor) equal for the normalised variables used in the model. The effect of the dimerisation process was also considered in the formulation of the equations, which leads to stoichiometric coefficients of value two in Equation (3.3). Following the ideas proposed in [14], a possible delay, *τ*, was included in the rate that describes the deactivation of *DpS*_{
nc
} in the nucleus and subsequent translocation to the cytosol.

The experimental data available were converted into the normalised scale [5] defined in Sup. Mat. A1 [see Additional file 1]. Data from two replicate experiments were available and used as independent experiments for parameter estimation. In these experiments, the cytoplasmic concentration of activated STAT5, [pSTAT5_{cyt}], and the concentration of activated EpoR/JAK2 complex, [pEpoR], were measured in several time points. Additionally, the extracellular concentration of Epo, [Epo], was measured in an independent equivalent experiment. Additional algebraic equations, reflecting the relation between the measured quantities (observables) and the variables, were defined in the model:

[Epo] = *Epo* [pEpoR] = *pEpJ* [pSTAT5_{cyt}] = 2·*DpS*

The variables on the left-hand side represent the observables, while the right-hand side represent the variables considered in the model. Further details on data processing are discused in Sup. Mat. A2 [see Additional file 1]. Several models were tested, including models with fixed predefined integer kinetic orders, with and without time delay in the STAT5 cycling, and several models with an increasing number of variable non-integer kinetic orders. In this latter case, variable non-integer kinetic orders related to terms representing a higher simplification of the dynamics were tested first (EpoR and STAT5 activation) and after that, the number of terms with non-integer kinetic order was systematically increased. The best compromise between an appropriate data fitting and a suitable number of parameters to be estimated was a model with fixed integer kinetic orders and no explicit time-delay for *DpS*_{
nc
} deactivation. Note, that although the chosen model here has no time delay as in [14], there is no contradiction between the two models. We do not represent the time delay implicitly but here we have a new state variable *DpS*_{
nc
} representing the fraction of activated STAT5 in the nucleus. The power-law model with non-integer kinetic orders for the term describing EpoR/JAK2 complex activation fitted the available data better (the objective function, *F*_{
obj
}, is 20% smaller than in the kinetic model), but the limited improvement obtained did not justify the choice of this model and we therefore followed the parsimony principle and selected the simpler and yet satisfactory model with fixed predefined integer kinetic orders. The procedure used and the different options analysed are discussed in Sup. Mat. A2 [see Additional file 1]. In all cases, parameters were computed for biologically relevant intervals in the parameter space using a genetic algorithm. The parameter values for the chosen model as well as the initial conditions used for parameter estimation are summarised in Table 1.

The model trajectories obtained for the chosen solution are depicted in Figure 2. The general behaviour of the system is reproduced despite the fact that the differences between the data obtained in the two replicates of the experiment are significant for some time points. We notice that the dynamics of phosphorylated EpoR is much better delimited by the experimental data and therefore the fit of the data is much clearer. The fit is also quantified with the objective function described in Equation 2 in Table A3 in Sup. Mat. A2 [see Additional file 1].

### Analysis of Amplification

To analyse how the signal is amplified through the pathway, we define the logarithmic amplification factor (*LA*) between two activated intermediates in a signalling pathway *X** and *Y** (with *Y** downstream in the pathway) with the following equation:

Where *T* is the duration of stimulation, *LA* is the logarithm of the ratio between the total productions of both intermediates during the signalling process considered. The total production of an intermediate is described by the integral of the net activation rate during the stimulation process. Considering this definition, a system amplifies between two steps in the pathway when *LA* is higher than zero. If *LA* is smaller than zero, the system provokes an attenuation of the signal. In the scale proposed, a value of one for *LA* implies that on average each molecule of *X** produces ten molecules of *Y**, while a value of minus one represents that ten molecules of *X** produce on average one molecule of *Y**. In Table 2 we propose a classification in significant intervals for the value of *LA* in terms of the ability of the system to amplify.

In order to analyse the responsiveness and the ability of the system to amplify signals, the performance of the system was analysed via mathematical simulation assuming three different conditions: sustained stimulation, transient stimulation and oscillatory stimulation by Epo. In case of a sustained stimulus, we analysed the response of the system in terms of the steady-state induced in the system for different values of constant Epo concentration in the extracellular medium, *Epo*_{
ss
}, from very low concentrations, *Epo*_{
ss
}= 10^{-8} units/ml, to concentrations up to tenfold the initial concentration used in the experiments performed, *Epo*_{
ss
}= 50.0 units/ml. Under normal conditions, the physiological serum concentration of Epo in mice is 7.9·10^{-3} units/ml [26]. We have computed the steady-state values of *DpS*_{
nc
}and they are shown in Figure 3. The logarithmic scale was used for the values of *Epo*_{
ss
}.

The figure shows clearly a sigmoidal behaviour in the log scale of *Epo*_{
ss
}. We can see that the system maximally responds to changes in concentration of Epo in the interval 10^{-3}-10^{-1}, which includes the normal range of Epo concentration in mouse serum (the behaviour is similar in case of *pEpJ* [see Additional file 1]). For smaller sustained stimuli the system is not significantly activated, while for intense stimuli the system reaches saturation and becomes virtually insensitive to any increase in the stimulus. Similar behaviour has been previously described as a typical feature of amplifying signal transduction pathways [7]. We analysed also the logarithmic amplification factor, *LA*, between the activated receptor, *pEpJ*, and the activated dimerised STAT5 in the nucleus, *DpS*_{
nc
}. The total production in steady-state is described by the integral of the net activation rate of the considered intermediate. If we use the suggested definition in this case, we obtain the following:

where ${V}_{pEpJ}^{+SS}$ is the steady-state value of ${V}_{pEpJ}^{+}$ and ${V}_{Dp{S}_{nc}}^{+SS}$ is the steady-state value of ${V}_{Dp{S}_{nc}}^{+}$ for the considered sustained concentration of *Epo*. *R*_{
STAT5/EJ
}is the ratio between the total amount of STAT5 and total amount of the EpoR at the cell surface. Since in our model we use normalised units, this ratio must be calculated in order to quantify the real value of amplification. The preliminary estimates for this ratio in the investigated cell line assigns a value around *R*_{
STAT5/EJ
}= 10 (data not shown). The inclusion of this factor is useful in order to quantify the numerical value of the amplification (it is acting as a translation factor in the *LA*-axis), but the rest of properties of the *LA* function (shape of the curve, scale, position of maxima and minima) depend on the intrinsic properties of the model and are not affected by the value of this ratio (*LA* is a logarithmic property; further explanations in Sup. Mat. A3 [see Additional file 1]). The definition in Equation (7.1) is valid for any kind of stimulus applied to the system (sustained, transient or oscillatory), while Equation (7.2) is the specific formulation for a sustained stimulus.

Figure 4 shows the results obtained for the considered interval of sustained stimulus. The logarithmic amplification factor has a value slightly higher than two (*LA* = 2) for all different values of sustained stimulation simulated, which implies an intense amplification of the signal. Thus, our model predicts that an activated receptor can on average activate and induce the nuclear translocation of up to 125 units of STAT5 before its deactivation. In addition, the maximal increase in amplification is in the interval 10^{-4}-10^{-1}; smaller stimuli produce a higher amplification, which implies that the system reacts sensitively to weak stimulation.

In case of a transient stimulus, the stimulation by Epo was characterised by two properties: the average duration of the stimulus, *T*_{
Epo
}, and the average value of the Epo concentration during transient stimulation, *Epo*_{
tr
}. In both cases we adapted the definitions proposed in [8] (a complete description is available in Sup. Mat. A4 [see Additional file 1]). Moreover, the variables considered for the analysis of a transient stimulation were the average amount of *pEpJ* (*pEpJ*_{
tr
}) and *DpS*_{
nc
}(*DpS*_{nc,tr}*)* during the stimulation, which are the transient stimulation "equivalent" for *Epo*_{
ss
}*and DpS*_{
ss
}. In the simulations, the initial values for the variables were consistent with the initial conditions (representing conditions after serum starvation) used in data fitting (Sup. Mat. A1 [see Additional file 1]). The stimuli analysed had a duration between 0.1 and 10^{4} minutes and an average stimulus concentration between 10^{-6} and 500 units/ml:

*T*_{
Epo
} ∈ [0.1, 10^{4}] min *Epo*_{
tr
}∈ [10^{-6},500] units/ml

Figure 5 shows the response of the system in *DpS*_{nc,tr}when transient stimulation of the system with the properties stated were simulated and analysed (Sup. Mat. A6 includes the figures on *pEpJ* [see Additional file 1]). The system reaches maximum average activation for stimuli with average intensity, *Epo*_{
tr
}, higher than 0.01 and average duration, *T*_{
Epo
}, between approximately 10 and 100 minutes. For longer stimulation, even at very high concentrations of Epo, there is a significant loss in *DpS*_{nc,tr}. For input signal longer than 500 minutes there is a partial recovery of the average signal intensity *DpS*_{nc,tr}and a plateau finally appears (which is even more apparent for *pEpJ*_{
tr
}[see Additional file 1]) for very long stimulation, which means that the system is not able to produce stronger responses for more intense stimuli. We suggest that this behaviour relates to the effect of limited recruitment of new EpoR to the plasma membrane, which is unable of faster recovery after very intense stimulation. We notice that is special features of the figure appear even for other computed solutions (data not shown) and seems to be a structural property of the model. In the range of the Epo pulse duration (*T*_{
Epo
}) between 1 and 100 minutes and for average stimulus concentration (*Epo*_{
tr
}) around the physiological values, [5.10^{-3}, 5.10^{-1}], the system shows maximum sensitivity to changes in the properties of the input signal. We furthermore investigated the responsiveness of the system predicted by the model with respect to the total amount of stimulus provided (*Epo*_{
tot
}), which in our simulations is represented by the product *Epo*_{
tot
}= *Epo*_{
tr
}· *T*_{
Epo
}(Sup. Mat. A4 contains further explanations [see Additional file 1]). Figure 6 shows our results, which condense the information depicted in Figure 5. The characteristic sigmoid-like behaviour of the model is also visualised in this figure. In addition, when a low stimulus is provided, the average response of the system (represented by *DpS*_{nc,tr}) does not depend on the features of the transient stimulus and stimulation with different duration and average intensity but the same total amount of Epo produces an almost identical response (Figure 6). In contrast, for a larger stimulus the response of the system will depend on the features of the input signal: different signals with the same total amount but different duration or average intensity will produce responses that can differ on average in ± 0.05 normalised units for *DpS*_{n,trc}, that is 50% of the average value (see solid black line in Figure 6). Interestingly, in case of medium and high *Epo*_{
tot
}, there is for any value of *DpS*_{
nc,tr
}an interval of values for *Epo*_{
tot
}in which signals that differ in all their properties will produce the same average output signal.

When we analysed the logarithmic amplification for a transient stimulus (Figure 7), we found that the range of values is again reduced but the average value is higher than two. The minimum amplification is reached for a transient stimulus with intermediate duration and high intensity, while the system shows a higher factor of amplification for short, weak stimulation in accordance with the results obtained for sustained stimulation. The maximum sensitivity of the amplification factor to changes in the signal is in the interval previously determined (*T*_{
Epo
}∈ [1,500 min], *Epo*_{
tr
}∈ [5.10^{-3}, 5.10^{-1}])

The dynamics of the system with oscillatory Epo concentration were also considered. For simulations, we did not choose perfect sinusoidal oscillatory signals but a design based on truncated triangular signals (Sup. Mat. A5 includes further explanations [see Additional file 1]). In [27] a physiological daily oscillation of the Epo concentration in the blood is described in which Epo is maintained at an almost stable value during daytime but reduces to half its value during the night. Two periods of transition are described that we adopted to describe the oscillatory signals used in our simulations. The oscillations in the input signal were characterised by two properties: the period of the oscillatory signal (T) which is the time between two consecutive maxima, and the average value of Epo during the oscillation (*Epo*_{
osc
}). The average value of *pEpJ* (*pEpJ*_{
osc
}) and *DpS*_{
nc
}(*DpS*_{nc,osc}) during the oscillation were defined and computed (Sup. Mat. A7 [see Additional file 1]). Our analysis suggested that the average of the signals for a number of periods between eight and twelve is sufficient to avoid the transient behaviour from the start of the simulations. We explored the performance of the system for oscillations with a period between half a minute and one day and with *Epo*_{
osc
}in the interval studied in the previous cases ([10^{-6},10] units/ml). Figure 8 shows the results obtained for *DpS*_{nc,osc}(Sup. Mat. A7 contains figures for *pEpJ*_{
osc
}[see Additional file 1]).

In Figure 8, we can see that the maximum response of the system for *DpS*_{
nc
} is reached for oscillatory signals with a period of 5–50 minutes and an average intensity between 0.1 and one (*DpS*_{
nc
} ≈ 0.15). This time interval coincides with the time scale suggested in [14] for the nucleocytoplasmic shuttling of STAT5, which could indicate a coupling between the frequency of the oscillations and the maximum performance of the pathway. For intense stimuli with longer period, the system reached a plateau in the average activated STAT5 in the nucleus at a value of half the maximum (*DpS*_{
nc
} ≈ 0.075). In case of the logarithmic amplification (Figure 9), a plateau with the highest value of amplification appears for oscillations with periods from short to long (T between 10–1000 min.) and weak stimuli (*Epo*_{
osc
} < 10^{–3}). Again, the maximum sensitivity to changes in the properties of the input signal appears for intermediate, physiologically feasible Epo concentration (*Epo*_{
osc
}) and for oscillations with a period T of 5–50 minutes. Intense stimulation with a short period produces a minimum in the amplification of the system, although the system is still strongly activated.

Interestingly, the properties of responsiveness and amplification predicted by the model with variable non-integer kinetic orders are very similar to the ones discussed here (Sup. Mat. A7 contains figures comparing both models [see Additional file 1]). The sigmoid-like curve on the responsiveness (*DpS*_{
nc
}) for sustained stimulation appears also for this model but is steeper and appears shifted one order of magnitude to higher values of sustained stimulus (Figure A7.1 [see Additional file 1]). Surprisingly, the logarithmic amplification (*LA*) in both models differs for very small sustained stimulus: the model with fixed predefined integer kinetic orders predicts a small variation on the values of the amplification through the whole interval of Epo values whereas the one with variable non-integer kinetic orders predicts an intense loss of amplification, which becomes negative (signal attenuation) for *Epo* ≤ 10^{-5} units/ml. Our analysis suggests that this feature is shared by the best power-law solutions (first ten of the selected power-law solutions, all them with steep sigmoidal curves). The other power-law solutions have a *LA* curve similar to the one obtained for the selected model (data not shown).

When transient stimulation is analysed, both models predict a similar shape for the surface accounting for the responsiveness and logarithmic gain of the system with respect to the features of the transient stimulation (*Epo*_{
tr
}, *T*_{
Epo
}), but again the sigmoidal curve for the model with variable non-integer kinetic orders is steeper and shifted (Sup. Mat Figure A7.2 [see Additional file 1]). The behaviour of the chosen model for oscillatory signals reproduces the features of the previously analysed signals (Sup. Mat Figure A7.3 [see Additional file 1]).

Finally, recent studies suggest the time-dependent recruitment of phosphatases like SHP-1 and other negative regulators like SOCS proteins as possible mechanisms to control cytokine responses, in particular in the pathway studied [28]. Any attempt to understand the behaviour of the JAK2-STAT5 pathway in depth requires to integrate the effect of these regulatory proteins in the mathematical model proposed and to generate the adequate experimental data to characterise these effects. The results that we show in this article should be confirmed or refined with this extended (feedback loop-controlled) model, additional experimental data and validation.

## Conclusion

A mathematical model based on ordinary differential equations and represented in power-law terms was developed for the JAK2-STAT5 pathway. Since the data available were rather limited, a strategy to reduce the complexity of the model was formulated following the strategy proposed in [15]. Using the available biological knowledge, several terms in the models were approximated by conventional kinetic terms with kinetic orders equal to one. With the remaining kinetic orders, a strategy of gradually increasing complexity in the structure of the model was used, which allowed a higher number of kinetic orders to be different to one for each iteration. The inclusion of a time delay in the STAT5 cycling was also considered. Although the best numerical fit to the data was obtained for a model with variable non-integer kinetic orders, a compromise between satisfactory reproducibility of experimental data and reduced complexity of the model lead to a simpler model with fixed predefined integer kinetic orders.

The responsiveness and amplification of the system were studied for sustained, transient, and oscillatory stimuli. Regarding responsiveness, we focussed on the value of activated STAT5 in the nucleus (*DpS*_{
nc
}), while the analysis of the fraction of activated receptors has been included in Sup. Mat. A7 [see Additional file 1]. To measure the ability of the system to amplify the signal we defined the logarithmic amplification (*LA*) as the logarithm of the ratio between the total amount of activated STAT5 in the nucleus (*DpS*_{
nc
}) and the activated receptor (*pEpJ*) during the analysed processes. A scale of values was set up in which positive values imply amplification of the signal, while negative values imply attenuation.

For the fraction of activated dimerised STAT5 in the nucleus (*DpS*_{
nc
}) and the logarithmic amplification (*LA*) the system seems especially sensitive in the range of physiological Epo concentrations [26]. Within this range changes in the properties of the stimulus for the three kinds of signals studied (sustained, transient, and oscillatory) provoke significant changes in the response of the system. For a stronger stimulus the system reaches saturation, while for a weaker stimulus there is no significant response (the system stays almost switched off). This could imply that the system is designed to have the maximum sensitivity to the intensity of the stimuli within the biologically feasible interval. The highest sensitivity of the system is not reached for transient stimulation or oscillatory input signals with very long duration (higher than one hundred minutes in both cases) but for intermediate values ([1,10^{2}] min). Thus, the system is set up to maximally respond to rapid changes in the environment but not to the long day and night oscillations registered under physiological conditions. Interestingly, the model predicts a small range of logarithmic amplification values (*LA≈* 2), which means that the average amount of STAT5 units activated per activated receptor remains very similar (around one hundred units) for a wide interval of Epo concentrations.

The system acts as a strong amplifier with respect to the scale defined in the three kinds of processes simulated and analysed. The model predicts that on average one activated dimeric EpoR can provoke the activation and subsequent nuclear translocation of approximately one hundred molecules of STAT5. Another interesting property is that the system seems to be more efficient when weak stimuli are applied.

For each model structure analysed, we stored a large collection of solutions and then analysed their properties. We found that, although the fitting to the data was not completely satisfactory in some of these solutions, a significant part of them showed similar properties related to the responsiveness and amplification of the system. Both type of models analysed (fixed predefined integer versus variable non-integer kinetic orders) show the characteristic sigmoid response curve, but there is a significant difference between them when the population of the solutions is considered. For the model with fixed integer kinetic orders there is a little variation between the 100 best solutions. On the other hand, for the the power-law model with variable non-integer kinetic orders there is a variety of solutions that reproduce sigmoid response curves including those of the model (as shown in Sup. Mat. A8 [see Additional file 1]). We think that this property relates to the results shown in [29], where the authors suggested that key properties of the biochemical networks, including signalling pathways, should be at least partially robust (in the sense of not-requiring "fine-tuning" of parameters) in order to ensure their proper functioning. With regard to this idea, we think that the sigmoid response curve obtained for the JAK2-STAT5 pathway mostly relates to the structure of the pathway and not to the precise values of the model parameters. In this case, the JAK2/STAT5 could be considered as a robust amplifier.

The final conclusion is that the JAK2/STAT5 system acts as an amplifier of the signal, which has the maximum sensitivity for input signals whose intensity coincides approximately with the physiological values, and reaches saturation for very intense and long stimulation. The general concepts, definitions and strategy of analysis proposed here could in principle be used to analyse the properties of any pathway once the dynamics of their regulatory proteins are measured.

## Methods

### Experimental techniques

Phoenix-eco cells were transiently transfected with the retroviral expression vector pMOWS-HA-EpoR by calcium-phosphate as previously described [30]. Six hours after transfection, Phoenix-eco cells were cultured in Iscove's Modified Dulbecco's Medium (IMDM) (Invitrogen) containing 50 μM β-mercaptoethanol (Invitrogen) and 30% fetal calf serum (Invitrogen). Virus-containing supernatant was harvested 24 h after transfection and filtered (0.45-μm filter, Corning). For spin infection, 10^{5} BaF3 cells were transduced with viral supernatants supplemented with 8 μg/ml Polybren (Sigma-Aldrich). BaF3 cells stably expressing HA-EpoR were selected in 1.5 μg/ml puromycin (Sigma-Aldrich) 48 hours after transduction.

For measuring activated Epo receptor and STAT5, HA-EpoR-BaF3 cells were starved in RPMI 1640 (Invitrogen) and 1 mg/ml BSA (Sigma-Aldrich) for 5 h. Cells were stimulated with 5 units/ml Epo (Janssen-Cilag) at 37°C for indicated times and 10^{7} cells per time point were lysed with NP-40 lysis buffer (1% NP-40, 150 mM NaCl, 20 mM Tris pH7.4, 10 mM NaF, 1 mM EDTA pH 8.0, 1 mM ZnCl_{2} pH4.0, 1 mM MgCl_{2}, 1 mM Na_{3}VO_{4}, 10 % glycerol) supplemented with aprotinin and AEBSF (Sigma-Aldrich). For immunoprecipitation, lysates were incubated with anti-EpoR antibodies (Santa Cruz) and anti-STAT5 antibodies (Santa Cruz). 40 ng of GST-EpoR and 36 ng GST-STAT5b was added as calibrator to each lysate as described [4]. Immunoprecipitated proteins were loaded in randomised fashion on SDS polyacrylamide gel as described [4], separated by electrophoresis and immunoblotted using anti-phosphotyrosine monoclonal antibody 4G10 (Upstate Biotechnology) and secondary HRP coupled anti-mouse antibody (Amersham Biosciences). Immunoblots were incubated with ECL substrate (Amersham Biosciences) for 1 min and exposed for 10 min on a LumiImager (Roche Diagnostics). LumiAnalyst software (Roche Diagnostics) was used for quantification. Antibodies were removed by treating the blots with β-mercaptoethanol (Sigma-Aldrich) and SDS (Serva) as described [31]. Reprobes were performed using anti-EpoR antibody (Santa Cruz) and anti-STAT5 antibody (Santa Cruz). Quantitative immunoblotting data was processed using GelInspector software [3]. Normalisers GST-EpoR for pEpoR and GST-STAT5b for pSTAT5 were used. For first estimates, csaps-splines were used with a smoothness of 0.3 for pEpoR and pSTAT5.

For measuring the extracellular Epo, an independent experiment was performed in which BaF3 cells stably expressing murine EpoR were starved in RPMI 1640 (Invitrogen) supplemented with 1 mg/ml BSA (Sigma-Aldrich) for 3 h and subsequently stimulated for up to 180 min with 5 units/ml [^{125}I]-Epo (Amersham Biosciences) at a density of 4 × 10^{6} cells/100 μl medium and 37°C. To separate free [^{125}I]-Epo from cell-associated [^{125}I]-Epo, cells were centrifuged through a layer of fetal calf serum (Invitrogen) and supernatants were measured in a gamma counter (Packard). Measurements were performed in triplicates. As control, cells were additionally incubated with an excess of unlabeled Epo (100 units/ml) (Janssen-Cilag), showing no decrease in free [^{125}I]-Epo in the medium.

### Mathematical modelling

The JAK2-STAT5 pathway has been modelled using a power-law representation [15]. As we discussed in Background, power-law models allow the representation of complex dynamics such as saturation-like behaviour, inhibition or cooperativity with simplified equations. Our strategy is to iteratively increase the complexity of our model by allowing variable non-integer kinetic orders for those processes where we do not have prior knowledge to determine the value of the kinetic order. We thus allow for a larger class of models, which includes the structure of conventional models based on mass-action kinetics (with fixed predefined integer kinetic orders) as a special case.

### Parameter estimation

In the present paper a genetic algorithm was used for parameter estimation. The algorithm has been adapted and optimised for power-law models. In the estimation process, each element of the population of solutions represents a point in the parameter value space. The initial population of solutions is generated through a random exploration of the search space, which is defined using feasible intervals of values for the variables (Sup. Mat. A2 [see Additional file 1]). The best individuals of the population are selected in the considered iteration based on the value of the following objective function [32]:

where *n*_{
exp
}is the number of experiments, *n*_{
var
}is the number of measured quantities (observables), *n*_{
tp
}the number of time points where each observable was measured, *X*_{
k, j
}(*t*_{
i
}) the value of the *j*^{th}observable at the *i*^{th}time point obtained after numerical integration of the solution for the *k*^{th}experiment. ${X}_{k,j}^{\mathrm{exp}}\left({t}_{i}\right)$ is the value of the *j*^{th}observable at the *i*^{th}time point measured in the *k*^{th}experiment. An additional fast-climbing stochastic algorithm is applied to the best solution each iteration of the algorithm. The stopping criterion is based on either a previously established maximum number of iterations or a minimum level of satisfaction for the objective function. Computations were performed on a Sun Fire V880 Server (four processors UltraSPARC-III, 1200 MHz, 8 MB cache each; RAM memory 32 GB). The algorithm for parameter estimation was implemented in Matlab R14 (The Mathworks, Inc. Natick, MA) running under SunOS 5.10.

## References

- 1.
Blume-Jensen P, Hunter T: Oncogenic kinase signalling. Nature. 2001, 411 (6835): 355-65. 10.1038/35077225

- 2.
Pawson T, Warner N: Oncogenic re-wiring of cellular signaling pathways. Oncogene. 2007, 26 (9): 1268-75. 10.1038/sj.onc.1210255

- 3.
Schilling M, Maiwald T, Bohl S, Kollmann M, Kreutz C, Timmer J, Klingmüller U: Computational processing and error reduction strategies for standardized quantitative data in biological networks. FEBS J. 2005, 272: 6400-6411. 10.1111/j.1742-4658.2005.05037.x

- 4.
Schilling M, Maiwald T, Bohl S, Kollmann M, Kreutz C, Timmer J, Klingmüller U: Quantitative data generation for systems biology: the impact of randomisation, calibrators and normalisers. Syst Biol. 2005, 152 (4): 193-200.

- 5.
Albeck JG, MacBeath G, White FM, Sorger PK, Lauffenburger DA, Gaudet S: Collecting and organizing systematic sets of protein data. Nat Rev Mol Cell Biol. 2006, 7 (11): 803-12. 10.1038/nrm2042

- 6.
Sourjik V, Berg HC: Receptor sensitivity in bacterial chemotaxis. Proc Natl Acad Sci. 2002, 9 (1): 123-7. 10.1073/pnas.011589998.

- 7.
Barkai N, Alon U, Leibler S: Robust amplification in adaptive signal transduction networks. CR Acad Sci Paris. 2001, 2 (IV): 1-7.

- 8.
Heinrich R, Neel BG, Rapoport TA: Mathematical models of protein kinase signal transduction. Mol Cell. 2002, 9 (5): 957-70. 10.1016/S1097-2765(02)00528-2

- 9.
Chaves M, Sontag ED, Dinerstein RJ: Optimal Length and Signal Amplification in Weakly Activated Signal Transduction Cascades. J Phys Chem B. 2004, 108: 15311-15320. 10.1021/jp048935f.

- 10.
Legewie S, Bluethgen N, Herzel H: Quantitative analysis of ultrasensitive responses. FEBS Journal. 2005, 272: 4071-4079. 10.1111/j.1742-4658.2005.04818.x

- 11.
Kisseleva T, Bhattacharya S, Braunstein J, Schindler CW: Signaling through the JAK/STAT pathway, recent advances and future challenges. Gene. 2002, 285: 1-24. 10.1016/S0378-1119(02)00398-0

- 12.
Aaronson DS, Horvath CM: A road map for those who don't know JAK-STAT. Science. 2002, 296 (5573): 1653-5. 10.1126/science.1071545

- 13.
Klingmüller U: The role of tyrosine phosphorylation in proliferation and maturation of erythroid progenitor cells–signals emanating from the erythropoietin receptor. Eur J Biochem. 1997, 249 (3): 637-47. 10.1111/j.1432-1033.1997.t01-1-00637.x

- 14.
Swameye I, Mueller TG, Timmer J, Sandra O, Klingmüller U: Identification of nucleocytoplasmic cycling as a remote sensor in cellular signaling by data-based modeling. Proc Natl Acad Sci. 2003, 100: 1028-1033. 10.1073/pnas.0237333100

- 15.
Vera J, Balsa-Canto E, Wellstead P, Banga JR, Wolkenhauer O: Power-Law Models of Signal Transduction Pathways. Cellular Signalling. 2007, 19: 1531-1541. 10.1016/j.cellsig.2007.01.029

- 16.
Kopelman R: Rate processes on fractals: theory, simulations and experiments. J Stat Phys. 1986, 42: 185-200. 10.1007/BF01010846.

- 17.
Kopelman R: Fractal Reaction Kinetics. Science. 1988, 241: 1620-1626. 10.1126/science.241.4873.1620

- 18.
Savageau MA: Development of fractal kinetic theory for enzyme-catalysed reactions and implications for the design of biochemical pathways. Biosystems. 1998, 47 (1–2): 9-36. 10.1016/S0303-2647(98)00020-3

- 19.
Brown N, Guoping L, Koszykowski ML: Mechanism reduction via principal component analysis. Int J Chem Kin. 1997, 29: 393-414. 10.1002/(SICI)1097-4601(1997)29:6<393::AID-KIN1>3.0.CO;2-P.

- 20.
Danø S, Madsen M, Schmidt H, Cedersund G: Reduction of a biochemical model with preservation of its basic dynamic properties. FEBS J. 2006, 273: 4862-4877. 10.1111/j.1742-4658.2006.05485.x

- 21.
Pohjanpalo H: System identifiability based on power series expansion of the solution. Mathematical Biosciences. 1978, 41: 21-33. 10.1016/0025-5564(78)90063-9.

- 22.
Jacquez JA, Greif P: Numerical parameter identifiability and estimability: Integrating identifiability, estimability, and optimal sampling design. Mathematical Biosciences. 1985, 77 (1–2): 201-227. 10.1016/0025-5564(85)90098-7.

- 23.
Vajda S, Rabitz H, Walter E, Lecourtier Y: Qualitative and quantitative identifiability analysis of nonlinear chemical kinetic models. Chem Eng Comm. 1989, 83: 191-219. 10.1080/00986448908940662.

- 24.
Hengl S, Kreutz C, Timmer J, Maiwald T: Data-based identifiability analysis of non-lineal dynamical models. Bioinformatics. 2007, 23 (19): 2612-2618. 10.1093/bioinformatics/btm382

- 25.
Gross AW, Lodish HF: Cellular trafficking and degradation of erythropoietin and novel erythropoiesis stimulating protein (NESP). J Biol Chem. 2006, 281 (4): 2024-32. 10.1074/jbc.M510493200

- 26.
Noe G, Riedel W, Kubanek B, Rich IN: An ELISA specific for murine erythropoietin. Brit J Haematol. 1999, 104: 838-840. 10.1046/j.1365-2141.1999.01273.x

- 27.
Jelkmann W: Molecular biology of erythropoietin. Intern Med. 2004, 43 (8): 649-59. 10.2169/internalmedicine.43.649

- 28.
Hilton DJ: Negative regulators of cytokine signal transduction. Cell Mol Life Sci. 1999, 55: 1568-77. 10.1007/s000180050396

- 29.
Barkai N, Leibler S: Robustness in simple biochemical networks. Nature. 1997, 387 (6636): 913-7. 10.1038/43199

- 30.
Ketteler R, Glaser S, Sandra O, Martens UM, Klingmüller U: Enhanced transgene expression in primitive hematopoietic progenitor cells and embryonic stem cells efficiently transduced by optimized retroviral hybrid vectors. Gene Ther. 2002, 9 (8): 477-87. 10.1038/sj.gt.3301653

- 31.
Klingmüller U, Lorenz U, Cantley LC, Neel BG, Lodish HF: Specific recruitment of SH-PTP1 to the erythropoietin receptor causes inactivation of JAK2 and termination of proliferative signals. Cell. 1995, 80 (5): 729-38. 10.1016/0092-8674(95)90351-8

- 32.
Wang FS, Ko CL, Voit EO: Kinetic Modeling Using S-systems and Lin-log Approaches. Biochemical Engineering Journal. 2006, 33: 238-247. 10.1016/j.bej.2006.11.002. 2007

## Acknowledgements

The authors acknowledge discussions with Christian Fleck and Thomas Maiwald from the University of Freiburg (Germany) and Edda Klipp from the Max-Planck Institute for Molecular Genetics (Berlin, Germany). This work was supported by the European Commission 6^{th} Framework program and as part of the COSBICS project under contract LSHG-CT-2004-512060 http://www.sbi.uni-rostock.de/cosbics. The contribution of Olaf Wolkenhauer was also supported by the German Federal Ministry of Education and Research (BMBF) grant 01GR0475 as part of the National Genome Research Network (NGFN-2) SMP-Protein.

## Author information

## Additional information

### Authors' contributions

JB, VB and ACP from the DKFZ (Heidelberg) carried out the experimental assays necessary to generate the data used in this work under the supervision of UK. JV from the University of Rostock designed the study, set up the mathematical model and performed the calculations concerning the analysis of responsiveness and amplification under the supervision of OW. JAH from University of La Laguna calibrated the model and generated the mathematical simulations used in this work under the supervision of NVTD. Finally all the authors including JT from University of Freiburg drafted the manuscript.

## Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

## Rights and permissions

## About this article

#### Received

#### Accepted

#### Published

#### DOI

### Keywords

- Kinetic Order
- Nucleocytoplasmic Shuttle
- Sustained Stimulation
- Transient Stimulation
- Transient Stimulus