Narrative-based computational modelling of the Gp130/JAK/STAT signalling pathway

Biological signalling pathways of even modest complexity cannot be comprehensively analysed within a feasible timescale by currently available experimental tools. However appropriate pathway models can be used to generate, explore and refine hypotheses guiding the formulation and prioritisation of experimental interventions. This has conventionally been approached by the use of models inspired by chemical kinetics and articulated mathematically in the form of ordinary differential equations. Recently an alternative approach, "molecules as computation", has been proposed in which a pathway is formulated as an executable computer programme [1, 2] which can be interrogated to determine the dynamic behaviour, robustness and parameter sensitivities of the model [3]. The outcomes of in silico experimentation on the computer model can then be used to inform the design of biological interventions in vitro. Among the various computational languages which have been used to model biochemical systems, we consider here process calculi (see for example [4–8], or [9] for a review of the process calculi approach). These calculi are formal languages which allow modellers to perform several kinds of analyses on the models (e.g. model-checking, causality and equivalence analysis).

One key challenge of this approach is the accurate description of biological pathways in the form of an executable computer language. From the biologist's perspective the formulation has to capture the biologically interesting features of the pathway and be readily understood by other biologists. From the computer scientist's perspective the formulation has to conform to the rules of formal methods in computer science: it must be logically precise and unambiguous. There is therefore a potential language gap between what the biologist understands and what the computer model encodes. We have recently described a high-level biologically-intuitive textual language in which the signalling pathway is articulated in the form of a narrative of events concerning the interactions between components located in different compartments [10]. This articulation of the pathway is then translated into an executable computer programme (written in the process calculus-based BlenX language [11]) for further analysis. In this paper we describe, develop and interrogate such an executable model of the gp130/JAK/STAT signalling pathway [12] using the Narrative Language approach and explore its predictions by in silico experimentation.

The gp130/JAK/STAT signalling pathway (see for example [12]) is the subject of significant clinical and biological interest, not least due to the key role it plays in human fertility, neuronal repair, haematological development and embryonic stem cell renewal [13]. Members of the gp130 cytokine family, such as LIF or OSM, bind to the common signal transducing receptor chain gp130 and a second signalling receptor LIFR or OSMR [14]. Homo- or hetero-dimerisation of gp130, LIFR and OSMR induces activation of the receptor associated kinase JAK which in turn phosphorylates the latent transcription factor STAT which, as a consequence, undergoes homo-dimerisation, translocates to the nucleus and activates the transcription of downstream gene targets (Figure 1). Several features of this pathway make it an attractive case study for a computer programme based modelling approach. A characteristic feature of the gp130/JAK/STAT system is the role of spatial confinement in which the transcription factor STAT undergoes nuclear/cytoplasmic shuttling which is regulated by JAK-mediated phosphorylation at the plasma membrane and T-Cell Protein Tyrosine Phosphatase (TC-PTP)-mediated de-phosphorylation in the nucleus [15, 16]. Aside from gp130/JAK activation by ligand the dynamics of the pathway can be regulated by a variety of mechanisms, which include STAT-mediated induction of Suppressor of Cytokine Signalling (SOCS) family proteins, which suppress JAK activation [17] and the Protein Inhibitor of Activated STAT (PIAS), an E3 family ubiquitin ligase which induces proteolytic degradation of phospho-STAT [18, 19]. Chronic nuclear STAT activation and/or JAK activation [20, 21] have been implicated as a predisposing event in a variety of tumour types indicating that pathway dynamics have significant impact on cell behaviour. Elucidating the relative influence of different pathway parameters on the activity of STAT will guide the evaluation of therapeutic interventions.

Much experimental data is available on the gp130/JAK/STAT pathway and several mathematical models have been developed based on ordinary differential equations (see for example [22–24]). We present a different approach to modelling this pathway based on formal computational methods. We use a novel textual language for modelling biochemical systems [10], and perform simulations through the BetaWB simulator [25], an existing tool based on Gillespie's stochastic simulation algorithm [26]. A key feature of the formal computational approach is the ability to rigorously explore, by in silico experimentation, the dynamic behaviour of the model to determine both the role of signal modifiers such as SOCS and the importance of parameter values. We validate our computational model by showing that it produces outputs which conform to those produced by experiment. We then perform in silico experiments on the model to determine first and second order parameter sensitivities and the effects of various types of pathway modulators. We show that the dynamic behaviour of the pathway is dominated by the rate of STAT de-phosphorylation and nuclear export and that these two variables result in bistable pathway behaviour when combined together in the model.

We used the N2BB tool to automatically translate the gp130/JAK/STAT pathway model into BlenX and we simulated the derived model using the BetaWB. Our intention in the simulations that follow was to firstly define the behaviour of the full model and then to study the behaviour of the model in response to perturbations such as modifying numerical values and/or removing components from the model. A particular feature of the narrative approach is that it is very simple to change parameters or modify the narrative of events in the process of model exploration.

Our aim is to define those features, which exert the dominating influence on the dynamic behaviour of the pathway in the model for evaluation by biological experiment. This is an example of using computer models to explore a wide range of scenarios in silico to guide the formulation of the more laborious and expensive laboratory-based experiments.

The time-evolution resulting from the simulation of the model is shown in the following pictures. Figure 2A reports the time-evolution of the full model, while Figure 2B–2D show the evolution when different inhibitors are removed (in silico genetics). The amounts of different STAT3 forms are plotted: monomeric cytoplasmic, dimeric cytoplasmic, monomeric nuclear, and dimeric nuclear.

Figure 2A shows that in the initial configuration STAT3 is present in monomeric form in the cytoplasm. As the system undergoes dynamic evolution, STAT3 is rapidly phosphorylated, dimerised and translocated into the nucleus. At the same time cytoplasmic STAT3 is dramatically reduced. The system reaches a plateau for some time and then slowly reverses as nuclear STAT3 levels fall and cytoplasmic levels rise. This behaviour conforms well to that obtained experimentally by observing the time-evolution of STAT3 phosphorylation (Additional file 9) and the relocation of STAT3 to the nuclear compartment in response to gp130 activation (Additional file 10). However it is important to note that, whilst the model and experimental data are in quantitative agreement for early time points (0–200 minutes), the experimental data reveals a faster rate of recovery of cytoplasmic STAT3 than predicted from the model over longer time periods. We were unable to accelerate the rate of relocalisation of STAT3 in the model by parameter variation (data not shown) which we interpret as indicating that the current model does not include biological processes that influence the rate of nuclear/cytoplasmic shuttling at later time points. In this context we remark that the current model takes no account of the induction of new gene expression [13] by gp130 signalling over this time period which could include components that influence the rate of nuclear/cytoplasmic transfer.

We conclude that the computer model derived from the original narrative is able to capture the dynamic behaviour of the real pathway, demonstrating the validity of the approach.

In the next phase we explored the dependency of the model on the presence of various components which were 'knocked out' by removal from the programme.

By comparing Figure 2A and 2B, in which we run the simulation in the absence of SOCS3, we observe that the effect of SOCS3 expression in response to STAT3 activation is to activate the slow attenuation process observed in the full model: removal of SOCS3 suppresses the delayed re-appearance of cytoplasmic STAT3 and the pathway exhibits prolonged and stable nuclear occupancy of STAT3.

The role PIAS3, which binds to phosphorylated nuclear STAT3 preventing it binding to DNA, is revealed by comparing Figure 2A and 2C. In this case removal of PIAS3 yields an initial increase of active STAT3, but the nuclear STAT3 signal attenuates with a faster time-course than in the full model, leading to an increase in cytoplasmic (and therefore inactive) STAT3.

Figure 2D shows the outcomes of removing the nuclear phosphatase TC-PTP from the model. In this case there is a rapid accumulation of nuclear STAT3 which reaches steady state and fails to attenuate. In this case cytoplasmic STAT3 is rapidly eliminated and does not re-appear.

Comparing the consequences of removing three different types of inhibitor from the model it is clear that each has characteristic temporal effects. The consequences of TC-PTP inhibition are significantly more rapid than removal of SOCS3 or PIAS3. Inhibition of SOCS3 and TC-PTP lead to prolonged and stable activated STAT3, whereas inhibition of PIAS3 accelerates the rate of activated STAT3 decay by inducing accumulation of cytoplasmic STAT3.

Single parameter sensitivity analysis

We next turned to analysis of the parameter sensitivities of the model. In this case we ran simulations in which individual parameters were systematically varied to observe the dependency of model behaviour on individual values. The aim here is to define those values, which have greatest impact on model behaviour.

We first examined the dependency of the model on the amount of STAT3 in the system. In this case we ran simulations containing different numbers of STAT3 molecules at initialisation (Figure 3). We observe that the duration of nuclear occupancy is relatively stable to STAT3 perturbation until it drops below a threshold value (in this case between 30–300 molecules) when it collapses into 'noise'. However the number of STAT3 molecules had a significant effect on peak amplitude where there was an approximately linear relationship between the maximum number of activated STAT3 molecules reached and the total number of molecules in the system.

We next examined the rate of TC-PTP de-phosphorylation, the amount of ligand, the rate of JAK kinase phosphorylation and the rate of nuclear export (Figure 4A–4D). The outcomes of these simulations show that the behaviour of the model is differentially dependent on particular values.

Slowing the rate of phospho-STAT3 de-phosphorylation exhibited significant effects on activated STAT3 amplitude and duration noted over the complete range analysed (Figure 4C). De-phosphorylation rates impacted on both the peak amplitude and duration of activated STAT3 indicating, as reported by others [32], that nuclear de-phosphorylation of activated STAT3 is an important determinant of signalling dynamics. Similar sensitivities were found on varying the rate of nuclear export (Figure 4D): although in this case nuclear export had no impact on peak signal amplitude and its main effect was on signal duration.

By contrast the model was relatively robust to variations in either the amount of ligand in the system or the rate of JAK activation (Figure 4A and 4B) where significant impacts on signal dynamics only become apparent at extreme values. Indeed, similar to the dependency on STAT3 numbers, it appears as though behaviour of the model is relatively robust to parameter changes in these processes over several orders of magnitude.

Multi-dimensional parameter sensitivity

The foregoing in silico experiments revealed that the dynamic behaviour of the gp130/JAK/STAT pathway is most sensitive to two parameters: the rate of nuclear STAT3 de-phosphorylation and the rate of STAT3 nuclear export. We were interested to learn if the model exhibited higher order dependencies when parameters were varied in combinations. Exhaustive implementation of this approach is currently computationally expensive. For this study we therefore chose to study the interaction between nuclear de-phosphorylation of STAT3 and export of de-phosphorylated STAT3. To this end 35 simulations were run in which each parameter was changed simultaneously (Figure 5). The results of this experiment were surprising. Instead of exhibiting a graded response across the whole parameter landscape, as might be predicted from the behaviour of each individual value, the system exhibits a bistable response in which, for the majority of conditions, the system is either constitutively activated or constitutively repressed. The system only exhibits sensitivity to parameter variation in a narrow region of values towards the middle of the ranges chosen. In this region both peak amplitude and signal duration were dependant upon the interplay between nuclear export and STAT3 de-phosphorylation.

Our purpose in this study was to explore the practical utility of the biological narrative approach for in silico exploration of a complex signalling pathway. We used the gp130/JAK/STAT pathway as a case study, and we showed that the biologist-specified narrative yields outcomes which conform well to experimental data. By 'experimenting' on the model by parameter exploration and component removal we were able to explore the influence of different elements on the dynamic behaviour of the pathway. Some of these experiments had already been performed in previous works [22–24], and the results we obtained show a good agreement with them; others provide novel insight on the pathway behaviour. Our analyses showed that the rate of STAT nuclear export and nuclear localised de-phosphorylation were key determinants of signalling dynamics. This conclusion is supported by in vitro experimental data in which inhibition of nuclear export by either drug treatment [33, 34] or mutation of the Nuclear Export Sequence [35] results in prolonged nuclear accumulation of phosphorylated STAT. The dependency of STAT import/export is also in agreement with existing mathematical models of the same pathway [22–24].

This outcome indicates that the model captures the dynamic behaviour of the real pathway well and encourages further exploration of the model into issues which would be resource intensive – or technically challenging – to address by biological experimentation. In particular we were interested in exploring the potential for interactions between parameters which are not currently addressable by biological experiment. We found, combining nuclear export and nuclear de-phosphorylation, that the two parameters interact strongly yielding a 'switch-like' behaviour. This type of modelling analysis may inform future considerations of multi-step mutagenesis or combination drug therapy scenarios in the gp130/JAK/STAT pathway.

The identified interaction between nuclear export and nuclear de-phosphorylation is novel, and the obtained behaviour was unpredicted; the insensitivity to the rate of JAK phosphorylation is another interesting outcome. These novel results reinforce our belief about the predictive power of this kind of modelling approach. In order to be validated, the hypotheses deriving from our simulation results should be tested by means of specific laboratory experiments. For instance, the characteristic dynamic shifts exhibited when the three inhibitors are removed from the system could be tested using appropriate mutant cells or by the use of specific pharmacological inhibitors. However, in general this might not be a simple task because it often involves the generation of such cells/drugs.

We have demonstrated in this study that a language can be used to describe biological signalling pathways in a way, which is formal and unambiguous for computational execution but more intuitive to the biologist compared with standard modelling languages. The approach exploits the particular advantages of the "molecule as computation" paradigm of Regev and Shapiro [1]: the resulting models are computable, relevant and understandable. There are also practical advantages to process calculus models in that they can be readily modified to explore different scenarios and interrogated using model-checking tools [3] to formally verify the model and explore its quantitative behaviour. Collectively the approach is therefore extensible in that, as new biological information on the pathway becomes available, it can be readily incorporated into the model.

In addition to the advantage of being simpler to use for modelling compared to directly specifying models in formal languages such as process calculi, the Narrative Language also provides modellers with a number of explicit features which capture some cardinal aspects of biological signalling pathways. It defines the temporal relationships between events (i.e. sequential, concurrent and competing events). It defines the location of proteins and the reaction volumes, and it deals with multi-compartmental models thereby making spatial location and confinement a central feature. Species in the model can exist in multiple states and locations. The Narrative Language can therefore be employed to model any biological process involving state transitions of different types; inter-molecular interactions; spatial confinement and temporal evolution. In order to assess the generalisability of the approach to model other biochemical systems, we have analysed a number of simple standard models using our and other existing simulation approaches and we have observed a good agreement between the results (results not reported here). Moreover, the biochemical entities, interactions and information which can be modelled in the Narrative Language are very similar to the ones modelled in the representation used in the NCI-Nature Pathway Interaction Database (PID) [36], a curated collection of biomolecular pathways represented in a graphical language. Considered the analogies between the PID and our proposal, we are currently developing a mapping between the PID representation and the Narrative Language. This translation would provide us with a significant number of well-understood pathways which can be directly simulated.

Finally the simulation results shown in this paper are obtained using N2BB, which translates Narrative Language models into BlenX models, but the translation of the Narrative Language into other current or future languages is possible. Thus the biological formulation of the pathway is separable from the computer method employed for simulation and analysis, allowing for interoperability with other tools for modelling and analysis of biochemical systems. For instance, a translation into the Bio-PEPA process algebra [37] is currently under development, and would allow modellers to describe models in the Narrative Language and study them by means of the various analysis methods usable for Bio-PEPA models (e.g. model-checking and mathematical techniques based on differential equations).