WO2002088383A2 - Method and system for identifying targets by nucleocytoplasmic cycling and use thereof - Google Patents

Abstract

The invention relates to a method, system and use thereof, in particular for identifying targets for efficient medical intervention and/or for predicting the effect of therapeutic agents and/or for dynamically modeling complex cell signaling networks and/or for remote sensing of cellular signaling. The method comprises the following steps: analyzing at least one cell signaling pathway, separating the cell signaling pathway into individual functional steps, transferring the individual steps into respective coupled differential equations comprising quantitative dynamical parameters by specifying the general differential equation (F) to describe the dynamics of individual steps of signal transduction such as covalent modifications, complex formation, cleavage, degradation, release and exchange.

Description

Method, System and Use thereof for Identifying Targets, for Efficient Medical Intervention and/or for Determining the Effect of Therapeutic Agents

The present invention relates to a method, a system and the use of the method and/or the system for identifying targets, for efficient medical intervention and/or for determining the effect of therapeutic agents and/or for remotely sensing cellular signaling by nucleocytoplasmic cycling.

Signaling pathways form complex intracellular networks that control proliferation, differentiation and survival. To understand these networks at the system level the dynamic interactions of individual components have to be examined which is facilitated by mathematical models (1 , 2). Previous attempts to model signaling pathways have been primarily based on qualitative data reflecting various interactions between the components and on simulations with ad hoc fixed parameters (3-6). However, to quantitatively predict the behavior of signaling pathways data driven models are required (7). To introduce meaningful simplifications and to establish a mechanism based model, detailed qualitative knowledge of a signaling pathway is necessary. A signaling pathway that has been studied in great detail is the JAK-STAT pathway (8-10). This pathway is involved in signaling through multiple cell surface receptors including hematopoietic cytokine receptors such as the erythropoietin receptor (EpoR) (11, 12). Signal transmission is initiated by ligand induced activation of a receptor bound Janus kinase (JAK) facilitating rapid tyrosine phosphorylation of the receptor cytoplasmic domain. The core module of the JAK-STAT pathway mediating rapid signal transduction from the cell surface receptor to the nucleus is represented by STAT-proteins that are tyrosine phosphorylated upon recruitment to the activated receptor and migrate as dimers to the nucleus where they stimulate the transcription of target genes (8-10). The molecular composition of signaling pathways has been studied in detail, but the dynamics of information processing is not understood. The dynamics of molecular signaling are, however, required for determining suitable targets, the activity of which must be modified for efficient medical treatment.

Many pharmaceutical compositions comprise drugs which act as modifiers of the activities of signalling molecules, such as ligands of receptor molecules. However, most conventional therapies are unspecific and could be replaced by more efficient and specific therapies. A critical issue for the efficacy and specificity of a therapy is to determine key targets which should be modified by a drug. If the key targets are known, a more effective drug can be identified by suitable screening assays or be designed based on known lead compounds. Suitable means for determining targets which are therapeutically valuable as described above have not been described yet but are nevertheless highly appreciated.

It is the object underlying the present invention to provide an improved method, an improved system and the use of the method and/or the system for identifying targets, for efficient medical intervention and/or for determining the effect of therapeutic agents and/or for remotely sensing cellular signaling by nucleocytoplasmic cycling.

This object is achieved with the subject-matter according to the claims.

Thanks to the present invention it has become possible to identify those targets, which are therapeutically valuable. Thus, the method of the present invention forms the basis for the development of new or improved drug therapies wherein the efficacy and specificity of the treatment is improved while the undesireable side effects are reduced or avoided. Advantageously, the method of the present invention can be integrated into a production process for said new or improved drugs. Preferably, said production process comprises at least the further step of producing, identifying and/or formulating a drug which efficiently modifies the activity of a target identified by the method of the invention in a therapeutically useful form. The present invention relates to a method to determine a mathematical model of the core module of the JAK-STAT signaling pathway based on experimental data. By in silico investigations of the fitted model the parameters of nuclear shuttling are determined as the most sensitive to perturbations and it is experimentally verified that inhibition of nuclear export results in a reduced transcriptional yield. The model reveals that a limited pool of STAT5 undergoes rapid nucleocytoplasmic cycles continuously coupling receptor activation and target gene transcription, thereby forming a remote sensor between nucleus and receptor. Thus, dynamic modeling of signaling pathways can promote functional understanding at the systems level.

Preferred but not limiting embodiments, examples and explanations are provided in the figures which show:

Fig.1 : Time course (points with error bars) and mathematical modeling (solid lines) of the STAT5 nucleocytoplasmic cycle. (A) Transient activation of EpoR and STAT5. Starved BaF3-EpoR cells (17) were left unstimulated (0) or were stimulated with 5 U/ml Epo for the time indicated. Cytoplasmic extracts were subjected to immunoprecipitation (IP) with anti-EpoR and anti-STAT5 antiserum, were analyzed by immunoblotting (IB) with the anti-phosphotyrosine (PTyr) antibody or reprobed with anti-STAT5 antiserum followed by enhanced chemiluminescence and detection using a Lumilmager. The Lumilmager files are displayed. (B) Linear interpolation of EpoR tyrosine phosphorylation as input function for the four dimensional differential equation. The time course of EpoR tyrosine phosphorylation in response to Epo stimulation was quantified with LumiAnalyst software and is displayed in arbitrary units. (C) and (D) show for the cytoplasmic tyrosine phosphorylated STAT5 and the total STAT5 pool in the cytoplasm the measured data in arbitrary units and the corresponding fit obtained with the linear model (C) and the model including nucleocytoplasmic cycling (D) (28).

Fig. 2: Predicting the behavior of the STAT5 nucleocytoplasmic cycle based on the dynamical parameters determined in the previous experiments. (A) Time course of EpoR tyrosine phosphorylation was used as input function to model (B) STAT5 tyrosine phosphorylation in the cytoplasm and (C) the total amount of cytoplasmic STAT5 in an independent experiment. Points with error bars indicate the experimental data whereas solid lines represent the mathematical modeling. The indicated error was determined based on duplicated measurements.

Fig. 3: In silico investigations. (A) Time courses of unobserved individual STAT5 populations. Depicted is the predicted quantitative behavior of unphosphorylated STAT5 (blue line), tyrosine phosphorylated STAT5 monomers (black line) and dimers (green line) in the cytoplasm and of cycling activated STAT5 molecules in the nucleus (red line). (B) Predicted the effect of parameter variations on target gene activation. As an indirect indicator for target gene activation the amount of nuclear activated STAT5 involved in cycling was determined by calculating the area under the red curve in (A). The effect of relative changes of the dynamical parameters k (black line), k₂ (blue line), 1 (green line), k_Λ (yellow line) and τ (red line) on the integrated area is shown.

Fig. 4: Effect of impaired nuclear export on the amount of activated STAT5 in the cytoplasm and the transcriptional yield. (A) Time courses of cytoplasmic STAT5 and EpoR phosphorylation as Lumilmager files (upper panels) and corresponding quantification of STAT5 phosphorylation. Starved BaF3-EpoR cells (16) untreated or pretreated for 30min with 10ng/ml leptomycin B were stimulated with 5 U/ml Epo for the time indicated. Cytoplasmic extracts were subjected to immunoprecipitation with anti-STAT5 antiserum followed by immunoblotting analysis with an anti- phosphotyrosine antibody. The delayed onset of signal activation in the presence of LMB is possibly due to the mild detergent function of LMB. Points with error bars represent experimental data quantified with LumiAnalyst software displayed in arbitrary units. The line indicates the model trajectories. (B) Reduced target gene production upon impaired nucleocytoplasmic cycling of STAT5. Starved BaF3-EpoR cells (10⁷ cells per time point) were left untreated or were pretreated with LMB and stimulated with 5 U/ml Epo for the time indicated. Cytoplasmic extracts were subjected to immunoprecipitation with anti-CIS antiserum raised against a GST-CIS fusion protein and analyzed by immunoblotting with anti-CIS antiserum followed by chemiluminescence detection. Quantification shown in the lower panel was performed using the LumiAnalyst software. The quantification of a representative experiment is shown in Boehringer Light Units (BLU). (C) Effect of nuclear export inhibition on activation of a STAT5 reporter gene. The STAT5 reporter construct pSac-CIS and as a control pSac-CIS-STAT lacking the STAT5 binding sites were introduced by electroporation into starved BaF3-EpoR cells. The cells were either left untreated or were pretreated with LMB and then unstimulated or were stimulated with 5 U/ml Epo. After 120min cell lysates were prepared and used for the determination of β^~-galactosidase activity. The results are displayed in relative light units (R.L.U.) and represent the mean of triplicate measurements ± SD normalized to protein content.

Fig. 5: Time course (points with error bars) and mathematical modeling (solid lines) of the STAT5 signaling pathway (a) fit with model 2, the feed-forward model with back-reaction (b) fit with model 4, the preferred model discussed above with additional back-reaction (c) fit with model 5, the preferred model discussed above with additional delay-distribution. The experimental data and the corresponding fits are shown for cytoplasmic tyrosine phosphorylated STAT5 (left panel) and total cytoplasmic STAT5 (right panel).

Preferred embodiments

To model the pathway, it is started with a parameterized model of the JAK-STAT core module reflecting the feed-forward information transfer from the cell surface to the nucleus and monitored ligand induced activation of STAT5 mediated by the EpoR. The individual steps of the core module were translated into four coupled differential equations describing the dynamics of the different STAT5 populations over time

(1) j , = -k_lx_lEpoR_A

(2) x₂ = -k₂x] + k_lx_lEpoR_A

) x₃ — k₃x₃ +—k₂x₂ (4) x₄ = +k₃x₃ ,

where the pool of cytoplasmic STAT5 is represented by *, (unphosphorylated STAT5), , (tyrosine phosphorylated monomeric STAT5), χ₃ (tyrosine phosphorylated dimeric STAT5) and x₄ is the nuclear STAT5. Since individual STAT5 populations are experimentally difficult to access it was measured in the cytoplasm the amount of tyrosine phosphorylated STAT5 (y_{ = k₅(x₂ +2x₃)) and the total amount of STAT5 (y₂ = k₆(x_l + x₂ +2x₃)). The scaling parameters k₅ and k₆ are introduced since only relative protein amounts are measured by quantitative experiments. As input function that determines the STAT5 response, Epo-induced tyrosine phosphorylation of the EpoR (EpoR_A , see Figs. 1A and B) was quantified (y₃ = k_ηEpoR_A). Tyrosine phosphorylation of the receptor reflects the extent of JAK2 activation and the effect of negative regulatory molecules including SOCS family members (13) and the tyrosine phosphatase SHP-1 (14). For simplification the protein inhibitor of activated STAT (PIAS) was not included in the current model. Potential cytoplasmic STAT-dephosphorylation is not considered, since the nucleus has been identified as the major compartment for STAT-dephosphorylation (15).

Due to the mathematical structure of the differential equations the scaling factor k_η can not be disentangled from the rate constant k_\ , and k₂ is coupled to , (0) (for detailed discussion of the identifiable parameter combinations see (16)). The identifiable parameter combinations were simultaneously estimated based on results from three independent experiments while the nuisance parameters k₅ ,k₆ and kη were determined separately for each experiment.

To estimate the parameters, time course experiments were performed with BaF3 cells expressing the wild type EpoR (BaF3-EpoR)(17), since these cells have been shown to efficiently activate the JAK-STAT pathway in response to Epo binding (18). Starved BaF3-EpoR cells were stimulated with 5 U/ml Epo for 2 to 60min. The amount of tyrosine phosphorylated EpoR and STAT5 present in the cytoplasmic lysates of these cells was determined by quantitative immunoblotting (see Fig. 1A) (19). The measurements shown in Fig. 1 revealed that tyrosine phosphorylation of the EpoR and STAT5 both increased rapidly upon Epo addition. Whereas after 10- 15min of Epo addition receptor tyrosine phosphorylation started to decline, STAT5 tyrosine phosphorylation remained at almost maximum levels until receptor phosphorylation dropped below threshold levels after 30min resulting in a plateau of STAT5 phosphorylation between 10-30 min. The total amount of STAT5 in the cytoplasm initially rapidly declined in response to Epo stimulation, but increased again after prolonged incubation times. Based on three independent experiments the parameters were estimated and the resulting time courses are displayed in Fig. 1C, showing that this model is not able to describe the experimental data.

Therefore, since it has been reported that dephosphorylated nuclear STAT- molecules can be relocated back into the cytoplasm (10, 15, 20, 21) the Eqs. (1 ,4) of the model are changed to include nuclear export (for detailed discussion of model derivation including tests for other models, see (16))

(1 ') x, = -k_lx_lEpoR_Λ + 2k₄x₃ ^T

* ) % ^{= —} Λ4 3 + K X ,

The superscript parameter τ represents the time STAT5 molecules reside in the nucleus. It is composed of the time STAT5 is bound to DNA and the time until dephosphorylated STAT5 leaves the nucleus (x₃ = x₃(t- τ)). To preserve mass conservation, the condition 1 ≥ k₄ has to hold. If the export rate k₄ is smaller than the import rate ^ , STAT5 accumulates in the nucleus. The time courses of the fitted model are displayed in Fig. 1 D. They describe all features of the experimental data including the plateau of phosphorylated STAT5 in the cytoplasm between 10 and 30 min. The parameters and their 1σ confidence intervals are: :, = 0.021 min^"1 (+0.004/- 0.003) £₇, k₂ = 2Λ6 min^"1mol^"1 (+1.7/- 1.0) / , (0), ^ = 0.1066 min^"1 (+0.03/- 0.022), k₄ = 0.10658 min^"1 (+0.0016/- 0.0024) and τ = 6A min (+0.5/- 2.6). As judged from the delay time τ , on average, a single STAT5 molecule remains in the nucleus for approximately 6min. To validate our mathematical model, quantitative data from an experiment showing atypical EpoR tyrosine phosphorylation was used as input function (Fig. 2A) to predict the development of cytoplasmic STAT5 in this experiment. Using the dynamical parameters k_x - k^ , τ and , (0) from the previous experiments and fitting only the nuisance parameters k₅ , k₆ and k_η, a curve was obtained that closely matched experimental data determined for both tyrosine phosphorylated and total STAT5 in the cytoplasm (Fig. 2B and C), thus demonstrating the predictive power of our mathematical model.

Based on the dynamic model the quantitative behavior of STAT5 populations can be estimated that are difficult to access experimentally (Fig. 3A). The simulation of the fitted model shows that upon activation of the system, tyrosine phosphorylated STAT5 monomers are formed, but rapidly converted into STAT5 dimers, providing an explanation for the technical difficulty of detecting monomeric tyrosine phosphorylated STAT5 experimentally. In the cytoplasm the amount of tyrosine phosphorylated STAT5 dimers reaches an initial maximum after approximately 7min and then declines due to nuclear translocation. Dephosphorylated STAT5 monomers are subsequently exported from the nucleus and rapidly re-phosphorylated at the activated receptor, thereby resulting in a second maximum of tyrosine phosphorylated STAT5 dimers in the cytoplasm after approximately 17 min. Therefore by mathematical modeling it can be revealed that the experimentally observed plateau of tyrosine phosphorylated STAT5 in the cytoplasm between 10 to 30min is the result of repetitive reactivation facilitated by nucleocytoplasmic cycling of STAT5. Fig. 3A shows that unphosphorylated STAT5 in the cytoplasm is limited and approaches zero approximately 9min after Epo stimulation. If, in a simulation, the nuclear delay time is increased to 30min, it is not possible to describe the plateau of tyrosine phosphorylated STAT5 in the cytoplasm, because STAT5 re-enters the pool of cytoplasmic STAT5 molecules when the EpoR is already in an inactive state (data not shown). Therefore, the overall nuclear sojourn time of activated STAT5 mediating target gene activation depends on efficient nuclear export and on the presence of activated EpoR. Data driven mathematical modeling identifies nucleocytoplasmic cycling as an essential behavior of the JAK-STAT core module and should be applicable to other STAT-family members since STAT-1 , a STAT5 related molecule, is dephosphorylated in the nucleus and then exported to the cytoplasm (15, 20, 21).

Based on the fitted model, the effect is evaluated of relative parameter changes on the amount of nuclear tyrosine phosphorylated STAT5 involved in cycling which should be proportional to target gene activation. This in silico investigation allows one to determine which step of the STAT5 nucleocytoplasmic cycle is most sensitive to perturbations. As indicated by the slope of the curves in Fig. 3B, changing the parameters of nuclear shuttling (^,Λ₄ and ? ) had the most significant effect on the amount of STAT5 in the nucleus whereas changes in the rate of monomer-dimer interconversion had a minor impact. By setting the nuclear delay time τ to infinite or alternatively k₄ = 0 , thus suppressing nuclear export, we could estimate that a single STAT-activation cycle contributes only 45 % of activated STAT5 in the nucleus compared to successive STAT nucleocytoplasmic cycles. In accordance with the prediction the experimental data presented in Fig. 4A shows that in the presence of the nuclear export inhibitor leptomycin B (LMB) the amount of tyrosine phosphorylated STAT5 in the cytoplasm declines more rapidly than in the absence of the inhibitor, whereas the time course of EpoR tyrosine phosphorylation was not reduced by the addition of LMB. Mathematical analysis of the experimental data revealed that the presence of LMB decreased the nuclear export rate (k₄) by 60 % (Fig. 4A). As judged from the relative parameter variation shown in Fig. 3B, this results within 60min in a 40 % reduction of activated STAT5 in the nucleus involved in cycling and thus should reduce the output of the system significantly. To experimentally test the prediction, the induction of the endogenous STAT5 target gene CIS (cytokine inducible SH2 domain containing protein) (22) was analyzed. As shown in Fig. 4A CIS was induced in BaF3-EpoR cells upon Epo stimulation and reached maximum levels after 60-90min. In agreement with the prediction the presence of LMB resulted in a 37% reduction of CIS after 60min of Epo addition and this increased to a 66% reduction at 90min. To confirm that this effect was mediated by a reduced efficiency of the JAK-STAT pathway, reporter gene assays (23) were performed (Fig. 4C) comparing the Epo induced induction of β-galactosidase by reporter gene vectors harboring the CIS promoter limited to the STAT5 binding sites (pSac-CIS) or harboring point mutations that inactivate the respective sites (pSac- CIS-STAT^") (24). In the absence of LMB the addition of Epo for 120min resulted in a 12-fold increase in β-galactosidase activity produced by the reporter gene vector pSac-CIS whereas Epo had no effect on the reporter gene vector pSac-CIS-STAT that lacks the STAT5 binding sites. The presence of LMB resulted in a general 2-fold reduction of β-galactosidase activity in the absence of Epo or STAT5 binding sites whereas the STAT5 specific yield of Epo induced β-galactosidase activity determined after 120min was reduced by 70%. Thus as predicted by our model impaired nuclear export of STAT5 results in a reduced yield of the system and the effect increases over time. Similarly, it was reported for STAT-1 that the presence of LMB or the use of a STAT1 mutant impaired in nuclear export caused a reduction in reporter gene activation (25). Therefore the dynamic model conclusively explains the counterintuitive observation that a partial block of nuclear export results in a reduced transcriptional yield.

It is therefore demonstrated that a quantitative data-specified model of a signaling pathway can be derived and applied to determine the dynamic behavior of the system in response to external and internal changes. The analysis of the JAK-STAT pathway has been focused on the uni-directional flow of information from the cell surface to the nucleus (8, 26). In contrast it is demonstrate by mathematical modeling that the yield of the system is determined by successive nucleocytoplasmic cycles. The insights in the dynamics of the pathway lead to the following insight at the systems level: By rapid recycling, STAT5 continuously senses receptor activation and couples the signal emanating from membrane receptors with target gene activation in the nucleus. This identifies nucleocytoplasmic cycling as a remote sensor that increases the possibilities for fine-tuned regulation. Therefore the functional analysis leads from a static description of signaling components to a dynamic understanding at the systems level. Investigation of additional modeling alternatives

To select an appropriate model describing the dynamics of the STAT5 signaling pathway it is necessary to go through a stepwise process of testing various models suggested by biochemical knowledge. We started with the assumption that the signaling pathway represents a feed forward cascade (model 1) delivering a signal from cell surface receptors to the nucleus, but as described in the manuscript this approach has to be discarded. To improve this model, we first included the possibility that STAT5 dimers posses a certain instability and disintegrate to the moπomeric form. Such back-reactions can occur, but it is a priori not clear to what extent they influence the dynamical behavior of this signaling pathway. Therefore, we investigated if a feed forward model with back-reaction (model 2) is able to describe the experimental data satisfactorily. In comparison to the simple feed forward model, the differential equations are extended by one additional term describing the back-reaction:

ι = -faxiEpoRA Vι = fea + 2a a)

Ϊ2 = +kiXlEpoRA - kz%2 + Zfc-jES 3 2 = fcs(-Gi + -G2 + 2£3) is = -føzs + OLS z - M$zs m = krEpoR

However, the fit of this model to the experimental data is similar to the fit of the simple cascade, see Supplemental Figure 1A, and therefore this model is also not able to describe the data.

Next we tested a model including the possibility that STATS is exported from the nucleus and reactivated in the cytoplasm. This model with nucleocytoplasmic cycling (model 3) leads to a satisfactory fit and is therefore described in detail in the manuscript, see Fig. ID.

To compare two models statistically, we compute twice the difference of the log-likelihood of both models which leads to likelihood ratio tests (1). When comparing model 2 with model 3 we obtain the test statistic L = 2(LR_$ — LRe>) = 830.2 which corresponds to a p-value of < 10^-5, where LRi refers to the log-likelihood of model i Hence, the null hypothesis that the dynamic parameter responsible for cycling fci is zero has to be rejected.

the distribution of the test statistic does not follow the conventional single χ^-distributioπ but a mixture of

since the estimated parameter under the alternative is on the boundary of the parameter space, i.e. may not be negative (2). The critical level of the 5% significance πiveau is 2.71.

Having established that a model with nucleocytoplasmic cycling is able to describe the dynamic behavior of the system, we investigated if the inclusion of the back-reaction into the model further improves the mathematical description. To test the improvement of the inclusion of the back-reaction in model 3, we extended our dynamical model accordingly, resulting in model 4. The new differential equations are now more complex since we introduce one additional parameter fc describing the back-reaction:

£ι = -hτiEpoRji - - 2kiX₃ t - r) y_ = A^(ϊ₂ + 2^₃)

£2 = +kιxιEpσRA -

+ £2 + Z3)

The resulting fit is displayed in Supplemental Figure IB. Again, to assess the significance of this new parameter, we calculate the test statistic LR — 2(LRι — LR3) = 0.70. This value corresponds to a p-value of 0.20 which indicates that this is not a significant improvement of the model.

Finally, we relaxed our assumption of a fixed sojourn time r for STAT5 in the nucleus. We investigated whether an introduction of a distribution of delays (model 5) improves the description of the experimental data significantly. The delay distribution was realized by folding the time course of z_$(t) with a Gaussian kernel with mean delay time TQ and window width w. The new time course of STAT5 leaving the nucleus x_$ t) is computed as

and the new differential equations are:

£1 = -kixiEpoRA + kix$ t) yi = fe(-«?2 + 2τs)

£2 = +kιXιEpσR —

V2 = ZI + Z₂ + 2£₃)

£3 = -fcs£3 + Q-δfø-cl ys = fy poRA

The resulting fit is displayed in Supplemental Figure 1C Using the same analysis as above we compute the test statistic as LR = 2(LR — LRs) = 0.55 corresponding to a p-value of 0.23. Thus we conclude that this model also implies an insignificant enlargement.

In summary, we conclude that the model with nucleocytoplasmic cycling but with neither the back-reaction nor the delay-distribution is the appropriate model. It represents a good trade off between the complexity of the mathematical description and its ability to represent the experimental data.

Parameter identifiability

Due to the observation function of the system which only allows limited measurements of the dynamical behavior, not all dynamical and nuisance parameters as well as the initial value xι(0) can be extracted from the data. This implies that the system is not fully identifiable (3). Therefore it is possible to obtain a dynamical system with less parameters consisting of the same observation function. The following transformation of variables leads to such an identical representation:

*s ⁿ = h ^r« ⁼ Therefore we are only able to extract single dynamical parameters fc&.fct and r and the parameter combinations h2Xi(Q),kι/k ,kQ/k2 and fo fø from the measured data. The new system of differential equations reads as follows: vι = -rιVχD -r 2nvs t - τ) yι = r&fø + 2v₃)

02 = +W .D - 0.5ϋ| Wt = re ( ?ι + v₂ + %) us = -rjus + vl ys = hrEpoR_Λ = D ύ₄ = - -rsϋ3 - r₄Vs(t - r).

It is noteworthy that the most interesting parameters, k^k± and r, which determine the cycling behavior are not affected by this transformation and may be extracted with help of the measured data and the πoπ-ideπtifiability of the parameters does not affect the model selection procedure described above.

Experimental Setup and Evaluation

Cellular signalling pathways, mediating receptor tivity to nuclear gene activation, are generally regarded as feed-forward cascades. We analyse measured data of a partially observed signalling pathway and address the question of possible feedback eyeing of involved biochemical components between nucleus and cytoplasm. First we address the question of cycling in general, starting from basic assumptions about the system. We reformulate the problem as a statistical test leading to likelihood ratio tests under nonsfcandard conditions. We find that the modelling approach without cycling is rejected. Afterwards, to differentiate two different transport mechanisms within the nucleus, we derive the appropriate dynamical models which lead to two systems of ordinary differential equations. These descri e the dynamical behaviour of all components of the pathway, lb compare both models we apply a statistical testing procedure based on bootstrap distributions. We find that one of both transport mechanisms leads to a dynamical model which is rejected while the other model is satisfactory.

Cells respond to their environment through complex signal transduction pathways that frequently begin at the cell membrane. In response to binding of specific ligands the receptor's cytoplasmic domain transmits the signal across the membrane and triggers a cascade of events in the cell. This often can be described by a series of biochemical reactions which may lead to gene expression. Current research on signalling pathways and metabolic networks is still mainly dedicated to answering questions concerning the qualitative behaviour of biological systems, e.g. showing the impact of certain biochemical components on a system or investigating interde endenties between variables (Bhalla and Iyen- gar, 1999; Fussenegger et al., 2000). While leading to a broader understanding in general, to fully comprehend these biochemical networks it is necessary to quantitatively investigate their dynamical behaviour (Koshland, 1998; Campbell, 1999; Zheng and Flavel, 2000; Endy and Brent, 2001; Downward, 2001). This inevitably leads to hypothesis testing to compare different models with statistical tests which are mostly applied under non-standard conditions in the finite sample case.

Starting with only a few assumptions about the underlying dynamics we first resolve the question of cycling, i.e. whether the biochemical reactant which activates the target genes may leave the nucleus to take part again in the signalling process. Afterwards we address the question of modelling tile cycling behaviour with statistical tests, comparing two different models which are derived from two different biochemical hypothesis about the exact transport mechanism. In both cases mathematical! models of the dynamical behaviour are tested for compatibility with the measured data.

The paper is organised as follows: In Section 2 we describe the biochemical background and formulate the mathematical properties based on a priori knowledge. Furthermore we describe the experimental setup focusing on questions related to partial observability of the system and on measurement noise. In Section 3 we reformulate the mathematical problem and resolve the question of cycling as a testing problem. Moreover we describe how to deal with s ecial nonstandard conditions and present the results from this investigation. In Section 4 we present two candidates of global parametric models which comply with the results from the previous section and test for compatibility with the data.

2 The biological model and the experimental setup

Signalling pathways play a dominant role in information processing within the cell. Upon binding of messenger molecules to cell surface receptors signal transmission is initiated within the cell leading to a cascade of biochemical reactions. In many cases this leads to migration of s ecific components into the nucleus where target genes are stimulated (Klingmuller et al., 1996). A signaling pathway that has been studied in great detail is the JAK-STAT pathway (Darnell, 1997; Pellegrini and Dusanter-Fourt, 1997) which can be activated by several receptors. Upon binding of erythropoietin (Epo) to the Epo receptor (EpoR) several biochemical reactions take place. In a first step the unphosphorylated STAT5 component (Signal Transducer and Activator of Transcription) is phosphorylated. In a next step, phosphorylated STAT5 molecules form dimers which are then able to enter the nucleus and to stimulate activation of target genes, see Fig.l (Has el et al., 1996). Hence, this signalling pathway can be described as a dynamical system consisting of one activation function, the time course of the activated EpoR, and four dynamical variables: unphosphorylated STAT5, phosphorylated STAT5, dimeriased STAT5 and STAT5 in the nucleus. Although it is known that STAT5 leaves the nucleus and can be used in future signalling cascades, the time scale of this release is unknown and hence the question of cycling of STAT5 during one signalling cascade arises (Haspel and Darnel!, 1999). This would mean that STAT5 leaving the nucleus can be phosphorylated anew by the Epo receptor, thus leading to a cycling behaviour of the STAT5 compo- nent.

Unfortunately, there is currently no measurement technique available which permits the direct measurement of any system variable. This leads to unobserved components which is a general problem in the analysis of biochemical systems. Only the activation function of EpoR, total amount of STAT5 in the cytoplasm and total amount of phosphorylated STATS can be measured over time. Even worse, all measurements are only proportional to the true amount introducing additional unknown nuisance scaling parameters. Additionally, measurements are corrupted by observational noise. Hence, we deal with partially observed noisy data. Since the measurement technique involves a counting process with many realizations, the error distribution can be approximated by a Gaussian distribution- In Fig.2 we display a typical measurement of all three observation functions. Error bars in Fug.2a are omitted since Epo receptor activation serves as an input function to the sy tem and is not modelled in a parametric way.

To increase the power of the test it is possible to use independent experiments. Then, the dynamical parameters governing the biochemical reactions of STATS are the same for every experiment while nuisance parameters differ.

To simplify the mathematical modelling we will use the following notation: E(t) is the time course of the activated EpoR at time t_s serving as input function to the system. xι(t)_tx?z (t) and -Cg(f-) are the amounts of unphosphorylated, phosphorylated and dimeri-εed STATS in the cytoplasm, xAt-) is the amount of STATS in the nucleus.

With a priori knowledge from biochemistry some properties of the system are known and we are able to derive mathematical descriptions of some biochemical reactions:

• The phosphorylation process of unphosphorylated STATS xι(t) which depends on the activated EpoR leads to the reduction of unphosphorylated STATS described by \pι ιE\. Release of STATS from the nucleus is modelled by [ t-G*,]. This is a first order approximation for many other transport mechanisms and therefore a rather general approach.

• The observation function of total amount of phosphorylated STATS in the cytoplasm is zι = s + atø).

• The observation function of total amount of STATS in the c toplasm is

• The observation function of activated Epo receptor is z₃ = pγE. The initial conditions are x₂(0) = a^(0) = ^(0) = 0, whereas a?ι(0) has to be estimated from the data. • Conservation of mass: = C = otmst.

• AE parameters are positive.

These properties are rather general and apply to many signalling pathways. Hence, the following analysis is applicable to other systems. It is important to note that in this first investigation we start from basic assumptions since we do not assume any knowledge about a mathematical model of the dimerization process (xz -+ xa ) or the transport of dimers into the nucleus (a?a → -&4).

3 Mathematical modelling and the formulation of the testing problem

In the following we will treat the problem only with one experiment but the same ine of thought can be applied to the multi-experiment case. Using a priori knowledge from Section 2, we obtain the following ordinary differential equations (ODE) and observation functions: i_t = —pι%ιE +-p4,x<ι = Ps(^2 + a^) (I)

*2 = +P1X1E - /fø) #2 = Ps (a?ι + aτ2 + a^ ) (2) i_s = +/(a?₂) - ff(a^) s = pτE =: D (3) , = +s^,(a^) ^— P-374, (4) with pi and H unknown dynamical parameter , a , pβ and pr unknown nuisance parameters, and f(xz) and gfas) unknown functions describing the biochemical reactions of dimerization and transport into the nucleus. Parameters pa and g will be later used to parameterise functions / and g. Using the mass conservation and the starting conditions we obtain:

_{Xl t) =} ^L- ^Α ₍₅₎

Ps Ps

- _f . it) (t) ,_βv

Pe Pa

Pa and differential equation (1) can be rewritten as:

Testing for cycling now amounts to the question whether parameter p* is different from zero. We interpret the measured time series as realizations of random variables and reformulate this question as a problem of model selection between two different linear regression models A and B denned as follows:

= + KiXu —

— K3X_4A +^■ K₄, (12) with

Xu = h(ti) X2i = Z2(h)E h) Xm = zι(k)E(k) Xu = *.(**) (14) and parameters

Kι i, K₃ = 4, Kt = ₄Psa?ι.(0). (15)

To discriminate models A and B, we apply the test statistic L = 2(L_B — LA) originating from the likelihood ratio test with L and L_B being defined as:

We assume all observations to be corrupted by white noise and obtain

The distribution of the test statistic L critically depends on the amount of information that can be drawn from the system. Before we deal with the actual setting of the recorded experimental data, we will investigate other realistic experimental settings and their influence on the distribution of the test statistic. First we assume X_ti to be realizations of independent random variables and only realizations of to be corrupted by observational noise. Models A and B originate from the same model class and are nested. Therefore, under standard conditions, we would assume L = 2(L_B — LA) to be distributed as a x² -distribution with Δtζf = άfs — #A degrees of freedom where Δcζf is the difference of degrees of freedom and άf_M denotes the degrees of freedom of model M. For the case with one experiment we get df = 2 since model B has additional parameters K₃ and

see Fig.3, curve A. In fact, due to the additional constraints for all parameters, K_? > 0, this distribution is different from the standard distribution (Self and Liang, 1987). The constraints lead to a distribution consisting of a mixture of χ² distributions so that

see Fig.3, curve B. The distribution incorporating the inequality constraints leads to a new 95% quantile qo. = 4.23. It is noteworthy to recall the 95% quantile under standard conditions which is £₉.35 = 5.99. Hence the constraints are important to consider to avoid a loss of power of the test.

Under more realistic conditions, random variables Xf are not independent since they originate from the ordinary differential equations (1-4). Moreover all observations are corrupted by measurement noise. To the best of our knowledge, it is then not possible to analytically describe the distribution of the test statistic L. Therefore we simulated this distribution with the number of data points being the same as in the experiment, the parameter values based on estimates with the smaller model and the error distribution taken from a priori biochemical knowledge. The result is displayed in Fig.3, curve C. The new conditions have a huge impact on the test statistic with the new 95% quantile being = 14.7.

In the next setting, we assume that only time series -£_&(£) and their temporal derivatives are directly accessible. Then, eqs. (11) may be used to compute the realizations of the X_?. This leads to a new distribution of the test statistic due to the nonlinear transformation Again, we simulate this distribution which displayed in Fig.3, curve D. We obtain a value of ^.95 = 19.0. Finally, we investigate the experimental setting under which the analysed data were recorded, i.e. observations are restricted to zι,z and 23 n. eqs. (1-3). In this case temporal derivatives are not accessible by measurement and have to be estimated from the data. Here, we use a method based on splines, Hanke and Scherzer (2001). In this case, random variables Y_tXχ_t . . . , t follow complex dependency structures and beforehand it is difficult to predict the impact on the test statistic. Again, we use a simulation to compute the distribution of the test statistic L which is displayed in Fig.3, curve E, and yields gϋ.₉₅ = 30.5.

When using the data set to calculate the value of the test statistic, we obtain a value of L = 78.7 which corresponds to a p-value of p < 10^-4. Therefore we conclude that model A is not sufficient to describe the data. In Fig .4 a and b, a typical fit for both models for the data set is shown. Model A is not able to fit the data while model B only misrepresents some time points in the beginning.

Model selection between models A and B may be also resolved by approximating the probability distribution of estimated parameter p,ι with a bootstrap procedure (Efron, 1982; Hinkley, 1984; Efron and Tibshirani, 1993). Drawing 10⁴ different bootstrap samples from the original data set and computing parameter estimates for every bootstrap sample, we obtain an approximation of the probability distribution for p as displayed in Fig.5, solid line. It is important to note that using this distribution to test if p<% is significantly greater than zero does produce false positive results. Our simulation reveals that the estimate of ₄ is biased, see Fig.5 dashed fine. This stems from the fact that it is not possible to obtain a maximum likelihood estimator due to the complex dependency structure of the random variables. Therefore, when applying bootstrap procedures in such cases, it is necessary to explicitly simulate the distribution of parameters estimates under the null hypothesis. Nevertheless, if we compare both distributions of i, under the null and with the actual data set, in our case we may conclude that parameter p,ι is significantly greater than the expected value under the null. This confirms our results from the previous analysis of the test statistic .

4 Dynamical modelling

To understand the dynamical behaviour of the signalling pathway in general, it is necessary to find an appropriate mathematical description for the temporal evolution of al dynamical variables. Hence, we look for a dynamical model in form of a deterministic differential equation. Such a model is helpful for gaining insight into the hidden mechanisms of the pathway and for designing future ex eriments since it can be used to predict the behaviour of the system under perturbations. Possible alterations can be analysed theoretically and numerically in a fast and cheap manner.

Naturally, all deterministic differential equations to model the dynamical behaviour are only approximations to truth. However, analysis of the system with a priori knowledge, e.g. from biochemistry, normally leaves only a few candidate models which can be compared with model selection procedures.

Due to the previous investigations the model class is already reduced to models with cyclic behaviour. In the following we use a priori biochemical knowledge to mathematically model the dimerization process and the transport into the nucleus. This will enable us later to focus on comparing different models of cycling. The dimerization process (scfe -+^• xs) is modelled as the reduction of phosphorylated STATS ) and production of dimer (373) by [pz-cf]- The transport into the nucleus (x₃ -+ x&) is described as the reduction of dimer (x ) and production of nuclear STATS (37 ) by \p₃x₃]. On the other hand there are two possible transport mechanisms governing the export to the cytoplasm (x<ι -+^■ xι) about which no reliable information is available:

(1) A compartmental approach: reduction of nucleitic STATS (^₄) and production of unphosphorylated STATS (x_%) is proportional to the amount of nucleitic STATS and can be modelled by j ^]-

(2) A model with delay: reduction of nucleitic STATS (x±) and production of unphosphorylated STATS ( ι) is exactly the amount of dimerized STATS (xs) that entered the nucleus at time ( — T). This mechanism can be modelled by [ aar-a — τ)\ and represents the idea that dimerized STATS binds to the DNA for a certain time, is dephosphorylated quickly, and may then enter the cytoplasm again.

It is important to note that the discrimination between these parametric models of the transport mechanisms (1) or (2) with help of the measured data offers the possibility to gain information about a mechanism which currently cannot be investigated directly by measurements. In the following, terms from model Mi (% = 1,2) will be denoted as [ ]« in the mathematical description and the differential equations and observation equations are:

xι = -piXiE [+ 43?4,]i [+pa&*(t - τ)_2 zι = s(a?2 + a?s) (19) = +PL IE -

«2 = Ps(a?i + la + i_a) (20) is = - asts + Pzx S3 = rE •= D (21)

*4 = + 3S3 [-Pt l]l [-P33?3(t - T)]2. (22)

Due to the observation function of the system which aEows only limited measurements of the dynamical behaviour, not all dynamical and nuisance parameters as well as starting value xι(0) can be extracted from the data Therefore it is not possible to identify all parameters of the system. This non-identifiability often occurs in dynamical models with unobserved components (Chappeil et al., 1990; Mutter et al., 2002). It is possible to obtain a dynamical system with less parameters consisting of the same observation function. Therefore it is necessary to apply a transformation of the differential equation to ensure that all transformed parameters are identifiable. Using the transformation

Vi = p₂a;ι n. = — r₃ = ps (23)

PT r₄ = p₄ r₅ = ^ r_s = ^, (24)

P2 P2 we obtain the foEowing differential equations: v_t = -r i D [+r₄Ui]ι [+r₃t^(έ - r)]₂ Vi = r_B (tfc + ^s) (25)

V2 = -rιVιD - υ$ yz = rg (t\ + v% + vg) (26) v^ = — r₃vs + ι^ γ$ = D (27) ϋ₄ = +r V₃ [-r₄v₄]ι [-rsvsit - τ)]z- (28)

Still, for all parameters r$ it must hold that r_$ > 0.

Maximum likelihood estimates of all unknowns are computed by minimising

where N is the number of time points and M is the number of observed components. y (ti) are the data points at time t-i and }/ ^{ ti,r,Vι(0)) is the solution of the ODE for model M with dynamical parameters r and starting value ι (0) . For details of the minimisation technique used see (Bock, 1981; Schittkowski, 1995; Timmer et al., 1998). When comparing models i and Mz to determine the mechanism of cycling, we deal with a model selection problem which is rather difficult since we compare non-nested nonlinear models in the finite sample case. Unfortunately, strict results only exist for the asymptotic case (Cox, 1962; Pesaran and Deaton, 1978; Vuong, 1989). Therefore, we use a bootstrap procedure to deal with this problem. In a first attempt we apply a procedure based on ( Wahrendorf et al., 1987; Hall and Wilson, 1991). We first compute the log likelihood of both mod- els and obtain the difference Lu = [Li - Lz\. Afterwards, we generate residual bootstrap samples Y* of the data Y and obtain the bootstrap distribution of S* = [£f₂ — ₁.₂J by computing the log likelihood for both models for every bootstrap sample. The bootstrap distribution S* now serves as an approximation of the distribution of [Lu — LQ] where LQ is the true value. Thus, it is possible to test the null hypothesis Ho that Mi [respectively ₂J is the true model and therefore LQ < 0 [respectivelyL₀ > 0]. Applying this procedure, we neither reject model i nor M₂- This is consistent to the results of (Hall and Wilson, 1991) who show that this testing procedure has little power. Therefore, to increase the power, we apply a procedure suggested by (Hal and Wilson, 1991). We first simulate separately under fitted model i to obtain residual bootstrap samples Y*¹. Afterwards, similar to the previous_^ procedure, we fit both models Mi and M and compute the likelihood ratio L\ = [Ll^l _n- LI¹]- Again we are able to obtain the bootstrap distribution oϊS^*x = [ ji - Lu] and as before S*¹ serves as an approximation of the distribution of [ ₁₂ - LQ]. The same procedure is done under fitted model M₂ leading to bootstrap samples Y*². Again we fit both models Mi and M_t and obtain the bootstrap distribution S*² = [Ll — Lu]. With these approximations we may compute signincance levels in the case of each statistics and we obtain a = 0.007 for model i and α = 0.26 for M₂- This indicates that Mz is satisfactory while i is not. In Fig. 6, we display the resulting fit for selected model M₂.

5 Discussion

We have shown the cycling behaviour of a signalling pathway with specific statistical tests. Since the methods we used are rather general, they provide a common framework for analysing such biochemical systems. In a first analysis, we investigated cycling in general with help of a model selection procedure in a linear regression setting. Due to inequality constraints for dynamical parameters and variables which are a common requirement for signalling pathways, likelihood ratio tests under nonstandard conditions are the appropriate technique when analysing these systems. We have shown that different experimental settings may have a huge impact on the distribution of the test statistic. Moreover in most cases it is necessary to simulate this distribution under the null hypotheses since analytical results are not available.

In a second investigation, we compared two dynamical models which consist of different cycling mechanisms. This lead to a model selection problem with non-nested nonlinear models in the finite sample case. A specialised bootstrap technique turned out to have enough power to discriminate both models. Due to the improvement of measurement techniques and the general importance of signalling pathways we expect the interest in modeling these systems to grow rapidly in the near future. Therefore it is necessary to provide the mathematical means for model selection procedures to avoid false positive results. Especially in the non-limit case with few data points statistical testing under nonstandard conditions remains a crucial problem.

Figs. 6 to 11 shall provide a better but not limiting understanding of the invention and/or preferred embodiments thereof:

Fig. 6.

Schematic representation of the dynamical behaviour of the signaling pathway. The unknown time scale of the cycling x^ -*^■ _t is investigated.

Fig. 7:

Measured time series of a typical data set: (a) amount of activated Epo receptor (b) amount of phosphorylated STATS (c) amount of total STATS. For visualisation purp ses, we connected the data points with lines.

Fig. 8:

Cumulative distributions of the test statistic for different settings of experimental conditions as described in the text. Included in the figure is the 0.95 significance level and the corresponding critical values.

Fig. 9:

Best fit for two different linear regression approaches and fg of eqs. (10) and (11) for data Y. (dotted line with points confers to /A respectively fs solid line with triangles confers to Y):

(a) the smaller model without cycling,

(b) the larger model with cycling.

Fig. 10: Approximated probabiity density for the estimated dynamical parameter p . Left curve: under the nul hypothesis, right curve: with a bootstrap procedure with 10⁴ bootstrap samples.

Fig. 11 :

Fit of the selected global model M for the experimental date. set.

References, the disclosure of which are all herewith incorporated by reference:

1. H. Kitano, MIT Press (2001).

2. H. Kitano, Science 295, 1662-4. (2002).

3. G. von Dassow, E. Meir, E. M. Munro, G. M. Odell, Nature 406, 188-92. (2000).

4. M. Fussenegger, J. E. Bailey, J. Varner, Nat Biotechnol 18, 768-74. (2000).

5. U. S. Bhalla, R. lyengar, Science 283, 381-7. (1999).

6. C. Y. Huang, J. E. Ferrell, Proc Natl Acad Sci U S A 9Z, 10078-83. (1996).

7. D. Endy, R. Brent, Nature 409, 391-395 (2001).

8. J. E. Darnell, Jr., Science 277, 1630-5 (1997).

9. C. M. Horvath, Trends Biochem Sci 25, 496-502. (2000).

10. K. Shuai, Oncogene 19, 2638-44. (2000).

11. S. N. Constantinescu, S. Ghaffari, H. F. Lodish, Trends Endochnol Metab 10, 18-23. (1999).

12. U. Klingmϋller, EurJ Biochem 249, 637-47. (1997).

13. R. Starr, D. J. Hilton, Bioessays 21 , 47-52. (1999).

14. U. Klingmϋller, U. Lorenz, L C. Cantley, B. G. Neel, H. F. Lodish, Cell 80, 729-38. (1995).

15. R. L Haspel, J. E. Darnell, Jr., Proc Natl Acad Sci U S A 96, 10188-93 (1999).

16. see Supplementary Material.

17. The wildtype EpoR was cloned into the retroviral expression vector pMOWS (27) and introduced into BaF3 cells by retroviral transduction. Cell lines stably expressing the EpoR (BaF3-EpoR) were selected in the presence of puromycin (18). Starved BaF3-EpoR cells were stimulated with 5 U/ml Epo (Janssen-Cilag) for the time indicated. For each time point 10⁷ cells were taken from the pool of cells and lysed by the addition of 2 x NP-40 lysis buffer to terminate the reaction. Cytoplasmic extracts were subjected to immunoprecipitation with anti-EpoR and anti-STAT5 antiserum (both Santa Cruz Biotechnology), were resolved on 15% SDS-PAGE and transferred to nitrocellulose membranes. The detection was performed by immunoblotting with the anti-phosphotyrosine (PTyr) antibody 4G10 (Upstate Biotechnology) or reprobed with anti-STAT5 antiserum followed by enhanced chemiluminescence and detection using a Lumilmager (Roche Diagnostics). For treatment with leptomycin B (LMB; a kind gift from Dr.M.Yoshida, University of Tokyo) starved BaF3-EpoR cells were incubated for 30min with 10 ng/ml LMB before Epo addition.

18. U. Klingmϋller, S. Bergelson, J. G. Hsiao, H. F. Lodish, Proc Natl Acad Sci U S A 93, 8324-8. (1996).

19. The immunoblots were performed under standardized conditions, were incubated with enhanced chemiluminescence substrate (Amersham) for 1min and were exposed for 10min on a Lumilmager (Roche Diagnostics). For quantifications the LumiAnalyst software was used applying the single band analysis package with automated lane and band identification, flat background correction and slant correction.

20. R. L. Haspel, M. Salditt-Georgieff, J. E. Darnell, Jr., Embo J 15, 6262-8 (1996).

21. K. M. McBride, C. McDonald, N. C. Reich, Embo J 19, 6196-206. (2000).

22. A. Yoshimura et al., EMBO J. 14, 2816-2826 (1995).

23. For transient transfection BaF3-EpoR cells were starved in RPMI 1640/1% FCS for 3 h and used for electroporation applying 2 pmol of b-galactosidase reporter plasmids pSac-CIS and pSac-CIS-STAT and 42 mg carrier plasmid. Electroporated cells were either left untreated or were pretreated with 10 ng/ml LMB and subsequently stimulated with 5 U/ml Epo. After 120min whole cell extracts were prepared and analyzed for b-galactosidase activity. The b-galactosidase activity determined by duplicate measurements was normalized to the protein content of the lysates determined by Bradford assay and is displayed in relative light units (R.L.U.).

24. The reporter plasmid pSac-CIS was generated by amplifying the CIS promoter from genomic DNA and inserting the Sacl/Xhol subfragment, comprising the two STAT5 binding sites and the authentic transcriptional start site, into pGL2 basic (Promega) and transferring the promoter cassette via Kpnl and Xhol into pCMVD (Clontech) using the EcoRI and Xhol restriction sites in the multiple cloning site. By overlap extension both consensus STAT5 binding sites were altered to AAxxxTTxxxxAAxxxTT resulting in pSac-CIS-STAT.

25. A. Begirt, T. Meyer, M. van Rossum, U. Vinkemeier, Proc Natl Acad Sci U S A 97, 10418-23 (2000).

26. S. Pellegrini, I. Dusanter-Fourt, EurJ Biochem 248, 615-33. (1997). 27. R. Ketteler, S. Glaser, O. Sandra, U. M. Martens, U. Klingmϋller, Gene Therapy, (2002), in press.

are the observed data for component j at time t_{ and yf (tι,xγ (t = 0),fc) is the model trajectory with starting value x_l(t = 0) and the parameters k ; note that x₂(t = o) = χ₃(t = o) = χ₄(t = o) = o. σ_ύ were estimated from duplicated measurements. Derivatives of χ² with respect to starting values and parameters were computed using an extrapolation scheme for finite differences as implemented by Ridder (30). Confidence intervals were computed with help of likelihood contours (31).

29. E. Hairer, S. P. Norsett, G. Wanner, in Springer series in Computational Mathematics. (Springer- Verlag, 1993).

30. C. J. F. Ridders, Advances in Engineering Software (1982), vol. 4.

31. W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. T. Flannery, Numerical recipes (Cambridge University Press, 1992).

References for the additional modeling alternatives:

1. D. Cox, C. Hinkley, Theoretical Statistics. Chapman & Hall, London, (1974).

2. S. G. Self, K. Y. Liang, . Am. Stat. Ass. 82, 605 (1987).

3. L. Ljung, System Identification: Theory For lite User. Prentice Hall, NJ. (1999).

References for the experimental setup and evaluation:

Bhalla, U. and lyengar, FL (1999). Emergent properties of networks of biological signalng pathways. Science, 283:381-387.

Bock, H. (1981). Numerical treatment of inverse problems in chemical reaction kinetics. In Ebert, K., Deuflhard, P., and J ger, W., editors, Modetlm of Chemical Reaction Systems, pages 102-125. Springer, Berlin.

Campbell, P. (1999). Can physics deliver another biological revolution? Nature, 397:89.

Chappell, M., Godfrey, K., and Vajda, S. (1990). Global identifiability of the parameters of nonlinear systems with specified inputs: a comparison of methods. Math. Biosci., 102:41-73.

Cox, D- R. (1962). Further results on tests of separate families of hypotheses. Journal of the Royd Statistical Society, Series B, 24:406-424.

Darnell, J. (1997). Stats and gene regulation. Science, 277:1630-1635.

Downward, I. (2001). The ins and outs of signalling. Nature, 411:759-762.

Efron, B. (1982). The Jack fe, the Bootstrap and Other Resampling Plans. Society for Industrial and Applied Mathematics.

Efron, B. and Tibshirani, FL (1993). An introduction to the bootstrap Chapman & Hall.

Endy, D. and Brent, R. (2001). Modelling cellular behaviour. Nature, 409:391- 395.

Fussenegger, M., Bailey, J., and J.A.Vamer (2000). A mathematical model of caspase function in apoptosis. Nature Biotechnol, 18:768-774.

Hall, P. and Wilson, S. (1991). Two guidelines for bootstrap hypothesis testing. Biometrics, 47:757-762.

Hanke, M. and Scherzer, . (2001). Inverse problems light: Numerical differentiation. Amer. Math. Monthly, (108): 512-^21.

Haspel, R. and Darnell, J. (1999). A nuclear protein phosphatase is required for the inactivation of STAT-1. Proc. Nat Acad. Science, 96:10188-10193.

Haspel, FL., Salditt-Georgieff, M., and Darnell, J. (1996). The rapid inactivation of nuclear tyrosine phosphorylated STAT1 depends upon a protein tyrosine phosphatase. Embo J, 15:6262-8.

Hinkley, D. (1984). Bootstrap methods. . Royd St t. Soc. B, 50:321-337. Klingmuller, U., Bergelson, S., Hsiao, ., and Lodish, H. (1996). Multiple tyrosine residues in the cytosolic domain of the eaythropoietin receptor promote activation of STATS. Proc Natl Acad Sci USA, 93:8324^8328.

Koshland, D. (1998). The era of pathway quantification. Science, 280:852.

Mϋller, T.G., Noykova, N., Gyllenberg, M., and T mmer, J. (2002). Parame- ter identification in dynamical models of anaerobic waste water treatment. accepted in Math. Biosciences.

Pellegrini, S. and Dusanter-Fourt, I. (1997). The structure, regulation and function of the janus kinase (J AK) and the signal transducers and activators of transcription (STATs). Eur. J. Biochem., 248:615-633.

Pesaran, M. H. and Deaton, A. S. (1978). Testing non-nested nonlinear regression models. Ecomπmtrica, 46(3) :677-694.

Schittkowski, K. (1995). Parameter Estimation in Differential Equations- In Agarwal, E., editor, Recent Trends in O tmwxatA Theory and Applications. WSSIAA- World Scientific Publishing Co., Singapore.

Self, S. G. and Liang, K. Y. (1987). Asymptotic properties of maximum Ekeli- hood estimators and likelihood ratio tests under nonstandard conditions. / Am Stat Ass, 82 (398): 605-610.

TΛmmer, J. , Mϋller, T. , and Melzer, W. (1998) . Numerical methods to determine calcium release fiux from calcium transients in muscle cells. Biophys. J., 74:1694-1707.

Vuong, Q. H. (1989). Likelihood ratio tests for model selection and non-nested hypotheses. Econσmetrica, 57 (2): 307-333.

Wahrendorf, J., Becher, H., and Brown, C. (1987). Bootstrap comparison of non-nested generalized linear models: applications in survival analysis and epidemiology. Appl. Statist, 36(1):72-81.

Zheng, T. and Flav , FL. (2000). Death by numbers. Nature Biotech., 18:717- 718.

Claims

Claims

1. Method, in particular for identifying targets for efficient medical intervention and/or for predicting the effect of therapeutic agents and/or for dynamically modeling complex cell signaling networks and/or for remote sensing of cellular signaling, with the following steps:

(a) analyzing at least one cell signaling pathway,

(b) separating the cell signaling pathway into individual functional steps,

(c) transferring the individual steps into respective coupled differential equations comprising quantitative dynamical parameters by specifying the general differential equation x = f(χ,p) to describe the dynamics of individual steps of signal transduction such as covalent modifications, complex formation, cleavage, degradation, release and/or exchange.

2. Method according to claim 1 , wherein it further comprises the step of specifying the dynamical parameters by quantitative time course determinations.

3. Method according to claim 1 or 2, wherein it further comprises the step of analyzing the sensitivity of parameter changes to identify targets for efficient medical intervention and/or predicting the effect of therapeutic agents by determining the respective parameter changes.

4. Method according to any one of the preceding claims, wherein the ligand induced JAK-STAT signaling pathway or constitutive activation of this pathway is analyzed in step (a).

5. Method according to any one of the preceding claims, wherein the Epo- induced JAK-STAT signaling pathway is analyzed in step (a).

6. Method according to claim 5, wherein in step (d) the dynamical parameters are determined by measuring ligand induced tyrosine phosphorylation of the EpoR and STAT-5 and the amount of total STAT-5 in the cytoplasm.

7. Method, in particular for identifying targets for efficient medical intervention and/or for predicting the effect of therapeutic agents and/or for dynamically modeling complex cell signaling networks, the method comprising the following steps:

(a) storing into a system coupled differential equations comprising quantitative dynamical parameters by specifying the general differential equation £ = /(3c, ) to describe the dynamics of individual steps of signal transduction such as covalent modifications, complex formation, cleavage, degradation, release and/or exchange, the equations corresponding to the dynamics of individual functional steps in a cell signaling pathway,

(b) determining the dynamical parameters,

(c) introducing the dynamical parameters into the system,

(d) computing the differential equations with the dynamical parameters and

(e) outputting the result of the computing.

8. Method according to claim 7, wherein in step (a) coupled differential equations corresponding to the dynamics of the ligand induced JAK-STAT signaling pathway and/or constitutive activation of this pathway and/or Epo-induced JAK-STAT signaling pathway are stored.

9. Method according to claim 8, wherein in step (a) the following coupled differential equations are stored:

(1) i, = -p_lx_iR_A + 2p₃x₃ ^T

(2) x₂ = -p₂x] + p_λx_xR_A

wherein the pool of cytoplasmic STAT-5 is represented by JC, (unphosphorylated STAT-5), x₂ (tyrosine phosphorylated monomeric STAT- 5), (tyrosine phosphorylated dimeric STAT-5) and x₄ is the shuttling nuclear STAT-5, wherein the superscript parameter τ represents the time STAT-5 molecules reside in the nucleus and R_A represents the ligand induced or constitutive tyrosine phosphorylation of a receptor or a kinase and wherein the dynamical parameters p_v p₂,p₃, p₄ and τ are independently determined.

10. Method according to claim 9, wherein for R_A EpoR_A is used as an input function for the Epo induced tyrosine phosphorylation of EpoR.

11. Method according to claim 10, wherein in step (b) of claim 9 the dynamical parameters p₁, p₂,p₃, p₄ and t are specified by quantitative time course determinations with BaF3 cells expressing the wild type EpoR, stimulating BaF3-EpoR with 5 U/ml Epo for 2 to 60 min and determining the amount of tyrosine phosphorylated EpoR and STAT-5 present in the cytoplasmic lysates of these cells by quantitative immunoblotting.

12. System, in particular for identifying targets for efficient medical intervention and/or for predicting the effect of therapeutic agents and/or for dynamically modeling complex cell signaling networks, the system being particularly obtained or obtainable by a method of any one of the preceding claims and comprising the following features:

(a) means for storing into a system coupled differential equations comprising quantitative dynamical parameters by specifying the general differential equations = f(χ,p) to describe the dynamics of individual steps of signal transduction such as covalent modifications, complex formation, cleavage, degradation, release and/or exchange, the equations corresponding to the dynamics of individual functional steps in a cell signaling pathway,

(b) means for determining the dynamical parameters,

(c) means for introducing the dynamical parameters into the system,

(d) means for computing the differential equations with the dynamical parameters and

(e) means for outputting the result of the computing means.

13. System according to claim 12, wherein the storing means store coupled differential equations corresponding to the dynamics of the ligand induced JAK-STAT signaling pathway and/or constitutive activation of this pathway and/or Epo-induced JAK-STAT signaling pathway.

14. System according to claim 13, wherein the storing means store the following coupled differential equations:

(1) x, = -p_xx_λR_A + 2p₃x₃

(2) x₂

+ p x R_A

(4) x₄ = -p₃x₃ ^τ + p₃x₃ ,

wherein the pool of cytoplasmic STAT-5 is represented by x, (unphosphorylated STAT-5), .^ (tyrosine phosphorylated monomeπ^'c STAT- 5), χ (tyrosine phosphorylated dimeric STAT-5) and x₄ is the shuttling nuclear STAT-5, wherein the superscript parameter τ represents the time STAT-5 molecules reside in the nucleus and R_A represents the ligand induced or constitutive tyrosine phosphorylation of a receptor or a kinase and wherein the dynamical parameters p_l, p₂,p₃, p₄ and τ are determined by independent experiments.

15. System according to claim 14, wherein for R_A EpoR_A is used as an input function for the Epo induced tyrosine phosphorylation of EpoR.

16. Computer program comprising program code means for performing the method of any one of claims 1 to 11 when the program is run on a computer.

17. Computer program product comprising program code means stored on a computer readable medium for performing the method of any one of claims 1 to 11 when said program product is run on a computer.

18. Computer system for performing the method of any one of claims 1 to 11.

19. Computer system according to claim 18, wherein the system comprises means for carrying out the method of any one of claims 1 to 11.

20. Use of a method according to any one of claims 1 to 11 for identifying targets for efficient medical intervention and/or for predicting the effect of therapeutic agents and/or for dynamically modeling complex cell signaling networks.

21. Use of a system according to any one of claims 12 to 14 for identifying targets for efficient medical intervention and/or for predicting the effect of therapeutic agents and/or for dynamically modeling complex cell signaling networks.

22. Targets obtained by a method, a system, a computer program, a computer program product, a computer system or a use according to any respective one of the preceding claims.

23. Method of producing a pharmaceutical composition comprising the steps of any one claims 1 to 11 and formulating the compound with a pharmaceutically acceptable carrier or diluent.