BACKGROUND
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  <cit>  which can be interrogated to determine the dynamic behaviour, robustness and parameter sensitivities of the model  <cit> .
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 .
these calculi are formal languages which allow modellers to perform several kinds of analyses on the models .
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  <cit> .
this articulation of the pathway is then translated into an executable computer programme  for further analysis.
in this paper we describe, develop and interrogate such an executable model of the gp130/jak/stat signalling pathway  <cit>  using the narrative language approach and explore its predictions by in silico experimentation.
the gp130/jak/stat signalling pathway  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  <cit> .
members of the gp <dig> cytokine family, such as lif or osm, bind to the common signal transducing receptor chain gp <dig> and a second signalling receptor lifr or osmr  <cit> .
homo- or hetero-dimerisation of gp <dig>  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 .
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 -mediated de-phosphorylation in the nucleus  <cit> .
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  family proteins, which suppress jak activation  <cit>  and the protein inhibitor of activated stat , an e <dig> family ubiquitin ligase which induces proteolytic degradation of phospho-stat  <cit> .
chronic nuclear stat activation and/or jak activation  <cit>  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 .
we present a different approach to modelling this pathway based on formal computational methods.
we use a novel textual language for modelling biochemical systems  <cit> , and perform simulations through the betawb simulator  <cit> , an existing tool based on gillespie's stochastic simulation algorithm  <cit> .
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.
methods
narrative modelling language
for modelling the gp130/jak/stat pathway, we use the narrative language, a high-level textual language which has been recently designed for modelling dynamically regulated and spatially-confined biochemical pathways.
here we introduce the main features of the narrative language; for a detailed description of the language and of its translation into an executable language, the reader is referred to  <cit>  and  <cit> .
the basic entities of the language are molecules  and sub-cellular locations .
in the language molecules can interact , undergo biochemical modification  and re-locate between compartments.
the time-dependent behaviour of the pathway is described in the form of a narrative of events involving these basic entities and functions, which imposes a temporal sequence and defines inter-dependencies and contingencies.
in the narrative approach each of the elements can denote 'real'  or abstract  entities.
in silico exploration of the pathway model is simply enabled by modifying the narrative description and/or changing parameter values.
a feature of this narrative approach is the use of reliability values associated with numerical parameters.
this is a percentage value which describes the reliability of the associated numerical value, and it can be used to distinguish between values that are highly certain because obtained from high quality biological experiment, and others which are inferred as a result of un-verified assumptions or 'guesswork' .
reliability values do not influence the behaviour of the programme but are annotations to inform use of the model.
in particular modellers can employ reliability values to identify parameter dependencies to be explored during model refinement.
model description and simulation
the narrative language model of the gp130/jak/stat signalling pathway is supplied as additional files  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig> and  <dig>  the molecular components we consider in the model are: two ligands , three membrane-bound receptors , one effector , and two inhibitors .
the receptor associated kinase jak and tc-ptp phosphatase are implicitly modelled.
additional file  <dig> shows the definition of the components involved in the pathway.
for each protein, the list of interaction/modification sites, states and locations are defined and initialised.
each receptor contains at least one ligand binding site , one binding site for socs <dig> inhibitor, and some phosphorylation sites.
moreover, receptors can be in dimeric state .
stat <dig> has one phosphorylation and four binding sites , and it can be monomeric or dimeric .
the initial quantities are also defined.
the number of ligands is calculated based on the known extracellular concentration .
four compartments are involved in the system : the exosol , the cell membrane , the cytosol , and the nucleus .
the compartment volumes are calculated based on the average cell radius and ratio between intra-cellular compartment volumes stated in  <cit> .
the number of spatial dimensions is used to distinguish between 2d compartments  and 3d ones.
additional file  <dig> shows the definition of the biochemical reactions occurring in the pathway.
reaction types  and rates  are defined here.
the reaction volume is the volume in which the reaction occurs , and it is used in deriving the actual rate at which the reaction occurs.
some of the reaction rates have been obtained from wet experiments, while others have been estimated based on information about similar reactions, or extracted from other models  <cit> ; reliability values are assigned to reaction rates and volumes.
finally, additional files  <dig>   <dig>   <dig>   <dig> and  <dig> show the definition of the narrative of events, which describes the evolution of the system; it is a sequence of basic events, each of which is a constrained textual description of a biochemical reaction involving at most two components.
moreover, events are grouped into processes for the sake of readability.
the defined events describe the binding/unbinding of ligand/receptor pairs, the downstream lif and osm pathways , the downstream stat <dig> pathway , and the inhibition mechanisms.
additional file  <dig> models the binding of ligands to receptors , and the inverse unbinding reaction .
for each event, the involved components and the occurring interaction are specified  together with the activating conditions , and the reference to the corresponding reaction in additional file  <dig>
additional file  <dig> models the dimerisation of pairs of receptor subunits to form receptor complexes , which is triggered by the binding of a ligand to one of the receptors  ; the dissociation of the receptor complexes is also modelled .
additional file  <dig> models the activation  of the receptor complexes, the binding of stat <dig> to a receptor complex, and the activation  of stat <dig>  once the receptor dimeric complex is formed, each receptor subunit  can phosphorylate on specific sites .
stat <dig> can bind on receptors' phosphorylated sites , and the binding of stat <dig> allows the phosphorylation of stat <dig> on site y <dig> .
additional file  <dig> models the unbinding of stat <dig> from receptor complexes, its homodimerisation, and nuclear/cytoplasmic shuttling .
once phosphorylated, stat <dig> can dissociate from the receptor complex ; the phosphorylated site allows stat <dig> to homodimerise .
when stat <dig> is in dimeric form, it can translocate into the nucleus  where it can carry out its specific functions : stat <dig> binds to the dna, activating the transcription of downstream gene targets.
nuclear stat <dig> is inactivated through de-phosphorylation by tc-ptp , which leads to its de-dimerisation , and its export to the cytoplasm , where stat <dig> can undergo additional cycles of activation.
additional file  <dig> models socs <dig> and pias <dig> inhibition mechanisms.
socs <dig> is produced by active stat <dig>  and degraded , and it acts by competing with stat <dig> in binding to receptors .
pias <dig> acts by binding to active nuclear stat <dig> .
we developed a tool  <cit> , n2bb, which implements an automatic translation of models described in the biologically-intuitive narrative language into executable computable models formulated in blenx  <cit> , a programming language inspired on the beta-binders process calculus  <cit> .
process calculi, originally developed for modelling mobile communicating systems, have recently been proposed as appropriate for simulating biological processes  <cit> , and they have proved themselves as powerful tools for dynamical modelling of complex biological systems  <cit> .
differently from differential equations, process calculi also allow for analysis of models .
the blenx model derived from the narrative language model is compatible with the betawb  <cit> , a collection of tools for modelling, simulating, and analysing blenx models.
hence, the model can be imported into the betawb designer, or directly simulated by means of the betawb simulator; the time-evolution of the simulation can be visualised by means of the betawb plotter or the snazer tool  <cit> .
for a detailed description of the blenx language and of the implementation of the simulator, see  <cit> .
cells, reagents and cytokines
mcf- <dig> human breast cancer cells were obtained from american type culture collection  and cultured as described  <cit> .
the human oncostatin m recombinant expression plasmid, pgex-3c-osm, was prepared, expressed and purified as described previously  <cit> .
western blot, immunofluorescence and data analysis
serum starved mcf- <dig> cells were stimulated with  <dig> ng/ml oncostatin m for increasing times  at 37°c.
for western blot analysis, cell lysates were prepared and analysed as described  <cit>  and monoclonal anti-phospho stat <dig>  and stat <dig>  antibodies used for immunodetection.
the density of the bands representing phospho-stat <dig> and stat <dig> were measured using imagej  <cit>  and expressed as the ratio of phospho-stat <dig> to stat <dig>  for immunofluorescence studies, mcf <dig> cells grown on coverslips were fixed with 4% paraformaldehyde , permeabilised with  <dig> % saponin solution  for  <dig> min and blocked for  <dig> hour in  <dig> % saponin solution  plus  <dig> % foetal calf serum.
cells were immunostained with monoclonal anti-stat <dig> antibody for  <dig> h at rt and incubated with texas red-conjugated secondary antibody containing hoechst  for  <dig> min at rt.
coverslips were mounted with  <dig> μl of mowiol solution  on the slide and observed under confocal microscope.
for localisation analysis, images captured were converted to greyscale and total stat <dig> fluorescence calculated from the sum of pixel density values using imagej.
nuclear stat <dig> fluorescence was calculated from selection of nuclear area  and cytoplasmic stat <dig> fluorescence calculated by subtracting nuclear staining from total cellular staining.
for each time-point analysis was performed on between  <dig> –  <dig> cells  and mean values ±  <dig> sd were calculated for total nuclear and cytoplasmic stat .
RESULTS
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 .
the amounts of different stat <dig> forms are plotted: monomeric cytoplasmic, dimeric cytoplasmic, monomeric nuclear, and dimeric nuclear.
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 socs <dig>  we observe that the effect of socs <dig> expression in response to stat <dig> activation is to activate the slow attenuation process observed in the full model: removal of socs <dig> suppresses the delayed re-appearance of cytoplasmic stat <dig> and the pathway exhibits prolonged and stable nuclear occupancy of stat <dig>
the role pias <dig>  which binds to phosphorylated nuclear stat <dig> preventing it binding to dna, is revealed by comparing figure 2a and 2c.
in this case removal of pias <dig> yields an initial increase of active stat <dig>  but the nuclear stat <dig> signal attenuates with a faster time-course than in the full model, leading to an increase in cytoplasmic  stat <dig>
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 socs <dig> or pias <dig>  inhibition of socs <dig> and tc-ptp lead to prolonged and stable activated stat <dig>  whereas inhibition of pias <dig> accelerates the rate of activated stat <dig> decay by inducing accumulation of cytoplasmic stat <dig>
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 stat <dig> in the system.
in this case we ran simulations containing different numbers of stat <dig> molecules at initialisation .
we observe that the duration of nuclear occupancy is relatively stable to stat <dig> perturbation until it drops below a threshold value  when it collapses into 'noise'.
however the number of stat <dig> molecules had a significant effect on peak amplitude where there was an approximately linear relationship between the maximum number of activated stat <dig> 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 .
the outcomes of these simulations show that the behaviour of the model is differentially dependent on particular values.
slowing the rate of phospho-stat <dig> de-phosphorylation exhibited significant effects on activated stat <dig> amplitude and duration noted over the complete range analysed .
de-phosphorylation rates impacted on both the peak amplitude and duration of activated stat <dig> indicating, as reported by others  <cit> , that nuclear de-phosphorylation of activated stat <dig> is an important determinant of signalling dynamics.
similar sensitivities were found on varying the rate of nuclear export : 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  where significant impacts on signal dynamics only become apparent at extreme values.
indeed, similar to the dependency on stat <dig> 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 stat <dig> de-phosphorylation and the rate of stat <dig> 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 stat <dig> and export of de-phosphorylated stat <dig>  to this end  <dig> simulations were run in which each parameter was changed simultaneously .
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 stat <dig> de-phosphorylation.
CONCLUSIONS
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  <cit> , 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  <cit>  or mutation of the nuclear export sequence  <cit>  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  <cit> .
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  <cit> : 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  <cit>  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 .
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 .
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   <cit> , 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  <cit>  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 .
authors' contributions
mlg designed the computational model, performed the simulation experiments, analysed the data and wrote the paper.
ad performed the laboratory experiments and analysed the data.
nud participated in designing the model, performed the laboratory experiments and analysed the data.
jkh participated in designing the model, analysed the data and wrote the paper.
cp participated in analysing the data and writing the paper.
all authors participated in conceiving the experiments, read and approved the final manuscript.
supplementary material
additional file 1
click here for file
additional file 2
click here for file
additional file 3
click here for file
additional file 4
click here for file
additional file 5
click here for file
additional file 6
click here for file
additional file 7
click here for file
additional file 8
click here for file
additional file 9
time-evolution of stat <dig> phosphorylation.
levels of phospho-stat <dig>  and total stat <dig>  after stimulation with oncostatin m.  densitometric analysis of immunoblots showing the time-course of phospho-stat3/stat <dig> ratios.
simulation time-course of phospho-stat3/stat <dig> ratios.
click here for file
additional file 10
nuclear localisation of stat <dig>   time-course of stat <dig> nuclear/cytoplasmic localisation.
images representative of stat <dig> localisation after stimulation with oncostatin m:  control;   <dig> minutes;   <dig> hours;   <dig> hours.
simulation time-course of stat <dig> nuclear/cytoplasmic localisation.
click here for file
acknowledgements
mlg is supported by the epsrc grant ep/e031439/ <dig> "stochastic process algebra for biochemical signalling pathway analysis".
ad is supported by the ec fp <dig> endotrack programme.
nud and jkh are supported by cancer research uk.
cp is supported by microsoft research.

 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 which can be interrogated to determine the dynamic behaviour, robustness and parameter sensitivities of the model. 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. these calculi are formal languages which allow modellers to perform several kinds of analyses on the models. 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. this articulation of the pathway is then translated into an executable computer programme for further analysis. in this paper we describe, develop and interrogate such an executable model of the gp130/jak/stat signalling pathway using the narrative language approach and explore its predictions by in silico experimentation. the gp130/jak/stat signalling pathway 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. members of the gp cytokine family, such as lif or osm, bind to the common signal transducing receptor chain gp and a second signalling receptor lifr or osmr. homo- or hetero-dimerisation of gp 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. 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 -mediated de-phosphorylation in the nucleus. 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 family proteins, which suppress jak activation and the protein inhibitor of activated stat, an e family ubiquitin ligase which induces proteolytic degradation of phospho-stat. chronic nuclear stat activation and/or jak activation 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. we present a different approach to modelling this pathway based on formal computational methods. we use a novel textual language for modelling biochemical systems, and perform simulations through the betawb simulator, an existing tool based on gillespie's stochastic simulation algorithm. 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. methods narrative modelling language for modelling the gp130/jak/stat pathway, we use the narrative language, a high-level textual language which has been recently designed for modelling dynamically regulated and spatially-confined biochemical pathways. here we introduce the main features of the narrative language; for a detailed description of the language and of its translation into an executable language, the reader is referred to and. the basic entities of the language are molecules and sub-cellular locations. in the language molecules can interact, undergo biochemical modification and re-locate between compartments. the time-dependent behaviour of the pathway is described in the form of a narrative of events involving these basic entities and functions, which imposes a temporal sequence and defines inter-dependencies and contingencies. in the narrative approach each of the elements can denote'real' or abstract entities. in silico exploration of the pathway model is simply enabled by modifying the narrative description and/or changing parameter values. a feature of this narrative approach is the use of reliability values associated with numerical parameters. this is a percentage value which describes the reliability of the associated numerical value, and it can be used to distinguish between values that are highly certain because obtained from high quality biological experiment, and others which are inferred as a result of un-verified assumptions or'guesswork'. reliability values do not influence the behaviour of the programme but are annotations to inform use of the model. in particular modellers can employ reliability values to identify parameter dependencies to be explored during model refinement. model description and simulation the narrative language model of the gp130/jak/stat signalling pathway is supplied as additional files and the molecular components we consider in the model are: two ligands, three membrane-bound receptors, one effector, and two inhibitors. the receptor associated kinase jak and tc-ptp phosphatase are implicitly modelled. additional file shows the definition of the components involved in the pathway. for each protein, the list of interaction/modification sites, states and locations are defined and initialised. each receptor contains at least one ligand binding site, one binding site for socs inhibitor, and some phosphorylation sites. moreover, receptors can be in dimeric state. stat has one phosphorylation and four binding sites, and it can be monomeric or dimeric. the initial quantities are also defined. the number of ligands is calculated based on the known extracellular concentration. four compartments are involved in the system: the exosol, the cell membrane, the cytosol, and the nucleus. the compartment volumes are calculated based on the average cell radius and ratio between intra-cellular compartment volumes stated in. the number of spatial dimensions is used to distinguish between 2d compartments and 3d ones. additional file shows the definition of the biochemical reactions occurring in the pathway. reaction types and rates are defined here. the reaction volume is the volume in which the reaction occurs, and it is used in deriving the actual rate at which the reaction occurs. some of the reaction rates have been obtained from wet experiments, while others have been estimated based on information about similar reactions, or extracted from other models; reliability values are assigned to reaction rates and volumes. finally, additional files and show the definition of the narrative of events, which describes the evolution of the system; it is a sequence of basic events, each of which is a constrained textual description of a biochemical reaction involving at most two components. moreover, events are grouped into processes for the sake of readability. the defined events describe the binding/unbinding of ligand/receptor pairs, the downstream lif and osm pathways, the downstream stat pathway, and the inhibition mechanisms. additional file models the binding of ligands to receptors, and the inverse unbinding reaction. for each event, the involved components and the occurring interaction are specified together with the activating conditions, and the reference to the corresponding reaction in additional file additional file models the dimerisation of pairs of receptor subunits to form receptor complexes, which is triggered by the binding of a ligand to one of the receptors; the dissociation of the receptor complexes is also modelled. additional file models the activation of the receptor complexes, the binding of stat to a receptor complex, and the activation of stat once the receptor dimeric complex is formed, each receptor subunit can phosphorylate on specific sites. stat can bind on receptors' phosphorylated sites, and the binding of stat allows the phosphorylation of stat on site y. additional file models the unbinding of stat from receptor complexes, its homodimerisation, and nuclear/cytoplasmic shuttling. once phosphorylated, stat can dissociate from the receptor complex; the phosphorylated site allows stat to homodimerise. when stat is in dimeric form, it can translocate into the nucleus where it can carry out its specific functions: stat binds to the dna, activating the transcription of downstream gene targets. nuclear stat is inactivated through de-phosphorylation by tc-ptp, which leads to its de-dimerisation, and its export to the cytoplasm, where stat can undergo additional cycles of activation. additional file models socs and pias inhibition mechanisms. socs is produced by active stat and degraded, and it acts by competing with stat in binding to receptors. pias acts by binding to active nuclear stat. we developed a tool, n2bb, which implements an automatic translation of models described in the biologically-intuitive narrative language into executable computable models formulated in blenx, a programming language inspired on the beta-binders process calculus. process calculi, originally developed for modelling mobile communicating systems, have recently been proposed as appropriate for simulating biological processes, and they have proved themselves as powerful tools for dynamical modelling of complex biological systems. differently from differential equations, process calculi also allow for analysis of models. the blenx model derived from the narrative language model is compatible with the betawb, a collection of tools for modelling, simulating, and analysing blenx models. hence, the model can be imported into the betawb designer, or directly simulated by means of the betawb simulator; the time-evolution of the simulation can be visualised by means of the betawb plotter or the snazer tool. for a detailed description of the blenx language and of the implementation of the simulator, see. cells, reagents and cytokines mcf- human breast cancer cells were obtained from american type culture collection and cultured as described. the human oncostatin m recombinant expression plasmid, pgex-3c-osm, was prepared, expressed and purified as described previously. western blot, immunofluorescence and data analysis serum starved mcf- cells were stimulated with ng/ml oncostatin m for increasing times at 37°c. for western blot analysis, cell lysates were prepared and analysed as described and monoclonal anti-phospho stat and stat antibodies used for immunodetection. the density of the bands representing phospho-stat and stat were measured using imagej and expressed as the ratio of phospho-stat to stat for immunofluorescence studies, mcf cells grown on coverslips were fixed with 4% paraformaldehyde, permeabilised with % saponin solution for min and blocked for hour in % saponin solution plus % foetal calf serum. cells were immunostained with monoclonal anti-stat antibody for h at rt and incubated with texas red-conjugated secondary antibody containing hoechst for min at rt. coverslips were mounted with μl of mowiol solution on the slide and observed under confocal microscope. for localisation analysis, images captured were converted to greyscale and total stat fluorescence calculated from the sum of pixel density values using imagej. nuclear stat fluorescence was calculated from selection of nuclear area and cytoplasmic stat fluorescence calculated by subtracting nuclear staining from total cellular staining. for each time-point analysis was performed on between – cells and mean values ± sd were calculated for total nuclear and cytoplasmic stat.  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. the amounts of different stat forms are plotted: monomeric cytoplasmic, dimeric cytoplasmic, monomeric nuclear, and dimeric nuclear. 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 socs we observe that the effect of socs expression in response to stat activation is to activate the slow attenuation process observed in the full model: removal of socs suppresses the delayed re-appearance of cytoplasmic stat and the pathway exhibits prolonged and stable nuclear occupancy of stat the role pias which binds to phosphorylated nuclear stat preventing it binding to dna, is revealed by comparing figure 2a and 2c. in this case removal of pias yields an initial increase of active stat but the nuclear stat signal attenuates with a faster time-course than in the full model, leading to an increase in cytoplasmic stat 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 socs or pias inhibition of socs and tc-ptp lead to prolonged and stable activated stat whereas inhibition of pias accelerates the rate of activated stat decay by inducing accumulation of cytoplasmic stat 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 stat in the system. in this case we ran simulations containing different numbers of stat molecules at initialisation. we observe that the duration of nuclear occupancy is relatively stable to stat perturbation until it drops below a threshold value when it collapses into'noise'. however the number of stat molecules had a significant effect on peak amplitude where there was an approximately linear relationship between the maximum number of activated stat 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. the outcomes of these simulations show that the behaviour of the model is differentially dependent on particular values. slowing the rate of phospho-stat de-phosphorylation exhibited significant effects on activated stat amplitude and duration noted over the complete range analysed. de-phosphorylation rates impacted on both the peak amplitude and duration of activated stat indicating, as reported by others, that nuclear de-phosphorylation of activated stat is an important determinant of signalling dynamics. similar sensitivities were found on varying the rate of nuclear export: 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 where significant impacts on signal dynamics only become apparent at extreme values. indeed, similar to the dependency on stat 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 stat de-phosphorylation and the rate of stat 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 stat and export of de-phosphorylated stat to this end simulations were run in which each parameter was changed simultaneously. 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 stat 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, 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 or mutation of the nuclear export sequence 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. 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: 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 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. 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. 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, 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 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. authors' contributions mlg designed the computational model, performed the simulation experiments, analysed the data and wrote the paper. ad performed the laboratory experiments and analysed the data. nud participated in designing the model, performed the laboratory experiments and analysed the data. jkh participated in designing the model, analysed the data and wrote the paper. cp participated in analysing the data and writing the paper. all authors participated in conceiving the experiments, read and approved the final manuscript. supplementary material additional file 1 click here for file additional file 2 click here for file additional file 3 click here for file additional file 4 click here for file additional file 5 click here for file additional file 6 click here for file additional file 7 click here for file additional file 8 click here for file additional file 9 time-evolution of stat phosphorylation. levels of phospho-stat and total stat after stimulation with oncostatin m. densitometric analysis of immunoblots showing the time-course of phospho-stat3/stat ratios. simulation time-course of phospho-stat3/stat ratios. click here for file additional file 10 nuclear localisation of stat time-course of stat nuclear/cytoplasmic localisation. images representative of stat localisation after stimulation with oncostatin m: control; minutes; hours; hours. simulation time-course of stat nuclear/cytoplasmic localisation. click here for file acknowledgements mlg is supported by the epsrc grant ep/e031439/ "stochastic process algebra for biochemical signalling pathway analysis". ad is supported by the ec fp endotrack programme. nud and jkh are supported by cancer research uk. cp is supported by microsoft research.
appropriately formulated quantitative computational models can support researchers in understanding the dynamic behaviour of biological pathways and support hypothesis formulation and selection by "in silico" experimentation. here we use a high-level narrative language in designing a computational model of the gp130/jak/stat signalling pathway and show that the model reproduces the dynamic behaviour of the pathway derived by biological observation. we then "experiment" on the model by simulation and sensitivity analysis to define those parameters which dominate the dynamic behaviour of the pathway. the results, which are shown to be in agreement with existing mathematical models of the gp130/jak/stat pathway, serve us as a form of validation of the model and of the approach itself. the described narrative model of the gp130/jak/stat pathway represents an interesting case study showing how, by using this approach, researchers can model biological systems without explicitly dealing with formal notations and mathematical expressions, nevertheless being able to obtain simulation and analysis results. we present the model and the sensitivity analysis results we have obtained, that allow us to identify the parameters which are most sensitive to perturbations. the model predicts that nuclear compartmentalisation and phosphorylation status of stat are key determinants of the pathway and that alternative mechanisms of signal attenuation exert their influence on different timescales.
 appropriately formulated quantitative computational models can support researchers in understanding the dynamic behaviour of biological pathways and support hypothesis formulation and selection by "in silico" experimentation. an obstacle to widespread adoption of this approach is the requirement to formulate a biological pathway as machine executable computer code. we have recently proposed a novel, biologically intuitive, narrative-style modelling language for biologists to formulate the pathway which is then automatically translated into an executable format and is, thus, usable for analysis via existing simulation techniques.  here we use a high-level narrative language in designing a computational model of the gp130/jak/stat signalling pathway and show that the model reproduces the dynamic behaviour of the pathway derived by biological observation. we then "experiment" on the model by simulation and sensitivity analysis to define those parameters which dominate the dynamic behaviour of the pathway. the model predicts that nuclear compartmentalisation and phosphorylation status of stat are key determinants of the pathway and that alternative mechanisms of signal attenuation exert their influence on different timescales.  the described narrative model of the gp130/jak/stat pathway represents an interesting case study showing how, by using this approach, researchers can model biological systems without explicitly dealing with formal notations and mathematical expressions, nevertheless being able to obtain simulation and analysis results. we present the model and the sensitivity analysis results we have obtained, that allow us to identify the parameters which are most sensitive to perturbations. the results, which are shown to be in agreement with existing mathematical models of the gp130/jak/stat pathway, serve us as a form of validation of the model and of the approach itself. 