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 JAK2STAT5 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 powerlaw 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 timecourse data. The final model, whether it is conventional model with fixed predefined integer kinetic orders or a model with variable noninteger 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 timeresolved quantitative data [35] 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 [710].
The Janus kinase – signal transducer and activator of transcription (JAKSTAT) pathways are one of the beststudied cell signalling pathways [11,12]. The JAK2STAT5 pathway is activated through various receptors, including the erythropoietin receptor (EpoR). Cytokineactivated 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).
Figure 1. Structure of the JAK2STAT5 pathway model proposed. Legend: Epo: concentration of Epo in the extracellular medium; EJ: fraction of nonactivated EpoR/JAK2 complex; pEpJ: fraction of activated EpoR/JAK2 complex; S: fraction of nonactivated and nondimerised STAT5 in the cytosol; DpS: fraction of activated and dimerised STAT5 in the cytosol; DpS_{nc}: fraction of activated and dimerised STAT5 inside the nucleus.
Previous work on databased mathematical modelling of the core module of the JAK2STAT5 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 powerlaw formulation that allows for noninteger 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 powerlaw models and conventional ODEs models used in systems biology is that kinetic orders can have noninteger values. There are two main reasons to allow noninteger kinetic orders: firstly, reactions in nonhomogenous environments lead to noninteger kinetic orders [1618] 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 powerlaw expressions, although other alternatives are possible [19,20]. In powerlaw 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 saturationlike 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 massaction kinetic model. By allowing noninteger, 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 powerlaw 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 nonparametric bootstrapbased 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 modelbased 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 JAK2STAT5 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 Epobound EpoR/JAK2 complex, pEpJ. The receptor activation was modelled with a term depending on the Epo concentration and the amount of nonactivated receptor complex at the plasma membrane 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 nonactivated EpoR/JAK2 complexes.
Three possible states have been considered for STAT5 in the model: i) nonactivated 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 (, 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.
Additional file 1. Supplementary material. The additional file contains further information about the model calibration process and some figures complementing the discussion and the model selection.
Format: PDF Size: 815KB Download file
This file can be viewed with: Adobe Acrobat Reader
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 powerlaw terms:
In this model we assume that, after the long period of starvation before starting the experiment, the system reaches a steadystate 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:
The variables on the lefthand side represent the observables, while the righthand 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 noninteger kinetic orders. In this latter case, variable noninteger 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 noninteger 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 timedelay 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 powerlaw model with noninteger 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.
Table 1. Values of the parameters in the selected solution
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].
Figure 2. Data fitting of the selected solution for the observables. A: activated EpoR, [pEpoR]. B: activated cytosolic STAT5, [pSTAT5_{cyt}]. Data from two replicates were used as independent experiments. The quantitative data obtained from the two experiments (data points represented as crosses for experiment 1 and points for experiment 2) are compared with the solution of the model identification process (lines). Experimental data were obtained from BaF3EpoR cells (proB cells exogenously expressing the murine EpoR) stimulated with 5 units/ml Epo and normalised as defined in Sup. Mat. A1 [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.
Table 2. Classification for the values of logarithmic amplification (LA) in significant intervals
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 steadystate 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 steadystate values of DpS_{nc }and they are shown in Figure 3. The logarithmic scale was used for the values of Epo_{ss}.
Figure 3. Steadystate values of DpS_{nc }(DpS_{nc,ss}) for different values of sustained stimulation on Epo (Epo_{ss}). The dashed black line indicates the physiological value for serum concentration of Epo (aprox. 7.9·10^{3 }units/ml), while the finely dashed line indicates the concentration of Epo used in the experiments performed (5 units/ml).
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 steadystate 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 is the steadystate value of and is the steadystate value of 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 LAaxis), 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.
Figure 4. Logarithmic amplification, LA (log. units, see Equation 7.1), for different values of sustained stimulation with Epo (Epo_{ss}). LA measures the signal amplification between pEpJ and DpS_{nc}. The dashed black line indicates the physiological value for serum concentration of Epo (aprox. 7.9·10^{3 }units/ml), while the finely dashed line indicates the concentration of Epo used in the experiments performed (5 units/ml).
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:
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 sigmoidlike 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.
Figure 5. Average fraction of dimerised phosphorylated STAT5 in the nucleus (DpS_{nc,tr}) (norm. units), during transient stimulation. The behaviour of the system was analysed for transient stimulation with an average duration of T_{Epo }∈ [0.1, 10^{4}] minutes, and an average concentration of Epo_{tr }∈ [10^{6},500] units/ml.
Figure 6. Average fraction of dimerised phosphorylated STAT5 in the nucleus (DpS_{nc,tr}) (norm. units versus the total amount of Epo (Epo_{tot}) provided during the transient stimulation. The system was analysed for transient stimulation with total amount of Epo Epo_{tot }∈ [10^{8},10^{6}] units/ml. The solid black line is used to highlight the different values of DpS_{nc,tr }obtained for a total amount of Epo of one hundred units. In this case the average of the system response ranges between 0.05 and 0.17 normalised units. The dashed black line is used to highlight how in the interval of intense stimulation input signals with totally different properties (Epo_{tot }∈ [10^{1},10^{5}] units/ml) produce output signal with identical average intensity DpS_{nc,tr}.
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}])
Figure 7. Logarithmic amplification, LA (log. units, see Equation 7.1), measured for different transient stimuli. LA measures the signal amplification between pEpJ and DpS_{nc}. The behaviour of the system was analysed for transient stimulation with an average duration of T_{Epo }∈ [0.1, 10^{4}] minutes, and an average concentration of Epo_{tr }∈ [10^{6},500] units/ml.
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]).
Figure 8. Average fraction of dimerised phosphorylated STAT5 in the nucleus (DpS_{nc,osc}) (norm. units) during oscillatory stimulation. The simulated data was averaged for several following periods. The behaviour of the system was analysed for oscillatory stimuli with an average period duration T ∈ [0.5, 1440] minutes, and an average concentration of Epo_{ocs }∈ [10^{6}, 10] units/ml. In all the simulations, the signal was averaged for 12 periods of the oscillatory stimulus.
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.
Figure 9. Logarithmic amplification, LA (log units, see Equation 7.1), measured for different oscillatory stimuli. LA measures the signal amplification between pEpJ and DpS_{nc}. The simulated data was averaged for several subsequent periods. The behaviour of the system was analysed for oscillatory stimuli with an average period duration T ∈ [0.5, 1440] minutes, and an average concentration of Epo_{osc }∈ [10^{6}, 10] units/ml. In all the simulations, the signal was averaged for 12 periods of the oscillatory stimulus.
Interestingly, the properties of responsiveness and amplification predicted by the model with variable noninteger kinetic orders are very similar to the ones discussed here (Sup. Mat. A7 contains figures comparing both models [see Additional file 1]). The sigmoidlike 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 noninteger 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 powerlaw solutions (first ten of the selected powerlaw solutions, all them with steep sigmoidal curves). The other powerlaw 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 noninteger 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 timedependent recruitment of phosphatases like SHP1 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 JAK2STAT5 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 loopcontrolled) model, additional experimental data and validation.
Conclusion
A mathematical model based on ordinary differential equations and represented in powerlaw terms was developed for the JAK2STAT5 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 noninteger 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 noninteger 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 powerlaw model with variable noninteger 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 notrequiring "finetuning" of parameters) in order to ensure their proper functioning. With regard to this idea, we think that the sigmoid response curve obtained for the JAK2STAT5 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
Phoenixeco cells were transiently transfected with the retroviral expression vector pMOWSHAEpoR by calciumphosphate as previously described [30]. Six hours after transfection, Phoenixeco cells were cultured in Iscove's Modified Dulbecco's Medium (IMDM) (Invitrogen) containing 50 μM βmercaptoethanol (Invitrogen) and 30% fetal calf serum (Invitrogen). Viruscontaining 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 (SigmaAldrich). BaF3 cells stably expressing HAEpoR were selected in 1.5 μg/ml puromycin (SigmaAldrich) 48 hours after transduction.
For measuring activated Epo receptor and STAT5, HAEpoRBaF3 cells were starved in RPMI 1640 (Invitrogen) and 1 mg/ml BSA (SigmaAldrich) for 5 h. Cells were stimulated with 5 units/ml Epo (JanssenCilag) at 37°C for indicated times and 10^{7 }cells per time point were lysed with NP40 lysis buffer (1% NP40, 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 (SigmaAldrich). For immunoprecipitation, lysates were incubated with antiEpoR antibodies (Santa Cruz) and antiSTAT5 antibodies (Santa Cruz). 40 ng of GSTEpoR and 36 ng GSTSTAT5b 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 antiphosphotyrosine monoclonal antibody 4G10 (Upstate Biotechnology) and secondary HRP coupled antimouse 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 (SigmaAldrich) and SDS (Serva) as described [31]. Reprobes were performed using antiEpoR antibody (Santa Cruz) and antiSTAT5 antibody (Santa Cruz). Quantitative immunoblotting data was processed using GelInspector software [3]. Normalisers GSTEpoR for pEpoR and GSTSTAT5b for pSTAT5 were used. For first estimates, csapssplines 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 (SigmaAldrich) 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 cellassociated [^{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) (JanssenCilag), showing no decrease in free [^{125}I]Epo in the medium.
Mathematical modelling
The JAK2STAT5 pathway has been modelled using a powerlaw representation [15]. As we discussed in Background, powerlaw models allow the representation of complex dynamics such as saturationlike behaviour, inhibition or cooperativity with simplified equations. Our strategy is to iteratively increase the complexity of our model by allowing variable noninteger 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 massaction 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 powerlaw 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. is the value of the j^{th }observable at the i^{th }time point measured in the k^{th }experiment. An additional fastclimbing 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 UltraSPARCIII, 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.
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.
Acknowledgements
The authors acknowledge discussions with Christian Fleck and Thomas Maiwald from the University of Freiburg (Germany) and Edda Klipp from the MaxPlanck 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 LSHGCT2004512060 http://www.sbi.unirostock.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 (NGFN2) SMPProtein.
References

BlumeJensen P, Hunter T: Oncogenic kinase signalling.
Nature 2001, 411(6835):35565. PubMed Abstract  Publisher Full Text

Pawson T, Warner N: Oncogenic rewiring of cellular signaling pathways.
Oncogene 2007, 26(9):126875. PubMed Abstract  Publisher Full Text

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:64006411. PubMed Abstract  Publisher Full Text

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):193200. PubMed Abstract

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):80312. PubMed Abstract  Publisher Full Text

Sourjik V, Berg HC: Receptor sensitivity in bacterial chemotaxis.
Proc Natl Acad Sci 2002, 9(1):1237. Publisher Full Text

Barkai N, Alon U, Leibler S: Robust amplification in adaptive signal transduction networks.

Heinrich R, Neel BG, Rapoport TA: Mathematical models of protein kinase signal transduction.
Mol Cell 2002, 9(5):95770. PubMed Abstract  Publisher Full Text

Chaves M, Sontag ED, Dinerstein RJ: Optimal Length and Signal Amplification in Weakly Activated Signal Transduction Cascades.
J Phys Chem B 2004, 108:1531115320. Publisher Full Text

Legewie S, Bluethgen N, Herzel H: Quantitative analysis of ultrasensitive responses.
FEBS Journal 2005, 272:40714079. PubMed Abstract  Publisher Full Text

Kisseleva T, Bhattacharya S, Braunstein J, Schindler CW: Signaling through the JAK/STAT pathway, recent advances and future challenges.
Gene 2002, 285:124. PubMed Abstract  Publisher Full Text

Aaronson DS, Horvath CM: A road map for those who don't know JAKSTAT.
Science 2002, 296(5573):16535. PubMed Abstract  Publisher Full Text

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):63747. PubMed Abstract  Publisher Full Text

Swameye I, Mueller TG, Timmer J, Sandra O, Klingmüller U: Identification of nucleocytoplasmic cycling as a remote sensor in cellular signaling by databased modeling.
Proc Natl Acad Sci 2003, 100:10281033. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Vera J, BalsaCanto E, Wellstead P, Banga JR, Wolkenhauer O: PowerLaw Models of Signal Transduction Pathways.
Cellular Signalling 2007, 19:15311541. PubMed Abstract  Publisher Full Text

Kopelman R: Rate processes on fractals: theory, simulations and experiments.
J Stat Phys 1986, 42:185200. Publisher Full Text

Kopelman R: Fractal Reaction Kinetics.
Science 1988, 241:16201626. PubMed Abstract  Publisher Full Text

Savageau MA: Development of fractal kinetic theory for enzymecatalysed reactions and implications for the design of biochemical pathways.
Biosystems 1998, 47(1–2):936. PubMed Abstract  Publisher Full Text

Brown N, Guoping L, Koszykowski ML: Mechanism reduction via principal component analysis.
Int J Chem Kin 1997, 29:393414. Publisher Full Text

Danø S, Madsen M, Schmidt H, Cedersund G: Reduction of a biochemical model with preservation of its basic dynamic properties.
FEBS J 2006, 273:48624877. PubMed Abstract  Publisher Full Text

Pohjanpalo H: System identifiability based on power series expansion of the solution.
Mathematical Biosciences 1978, 41:2133. Publisher Full Text

Jacquez JA, Greif P: Numerical parameter identifiability and estimability: Integrating identifiability, estimability, and optimal sampling design.
Mathematical Biosciences 1985, 77(1–2):201227. Publisher Full Text

Vajda S, Rabitz H, Walter E, Lecourtier Y: Qualitative and quantitative identifiability analysis of nonlinear chemical kinetic models.
Chem Eng Comm 1989, 83:191219. Publisher Full Text

Hengl S, Kreutz C, Timmer J, Maiwald T: Databased identifiability analysis of nonlineal dynamical models.
Bioinformatics 2007, 23(19):26122618. PubMed Abstract  Publisher Full Text

Gross AW, Lodish HF: Cellular trafficking and degradation of erythropoietin and novel erythropoiesis stimulating protein (NESP).
J Biol Chem 2006, 281(4):202432. PubMed Abstract  Publisher Full Text

Noe G, Riedel W, Kubanek B, Rich IN: An ELISA specific for murine erythropoietin.
Brit J Haematol 1999, 104:838840. PubMed Abstract  Publisher Full Text

Jelkmann W: Molecular biology of erythropoietin.
Intern Med 2004, 43(8):64959. PubMed Abstract  Publisher Full Text

Hilton DJ: Negative regulators of cytokine signal transduction.
Cell Mol Life Sci 1999, 55:156877. PubMed Abstract  Publisher Full Text

Barkai N, Leibler S: Robustness in simple biochemical networks.
Nature 1997, 387(6636):9137. PubMed Abstract  Publisher Full Text

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):47787. PubMed Abstract  Publisher Full Text

Klingmüller U, Lorenz U, Cantley LC, Neel BG, Lodish HF: Specific recruitment of SHPTP1 to the erythropoietin receptor causes inactivation of JAK2 and termination of proliferative signals.
Cell 1995, 80(5):72938. PubMed Abstract  Publisher Full Text

Wang FS, Ko CL, Voit EO: Kinetic Modeling Using Ssystems and Linlog Approaches.
Biochemical Engineering Journal 2006, 33:238247.
2007
Publisher Full Text