Stronger Computational Modelling of Signalling
Pathways Using Both Continuous and
Discrete-State Methods
Muffy Calder
1
, Adam Duguid
2
, Stephen Gilmore
2
, and Jane Hillston
2
1
Department of Computing Science, University of Glasgow, Glasgow G12 8QQ,
Scotland
2
Laboratory for Foundations of Computer Science, The University of Edinburgh,
Edinburgh EH9 3JZ, Scotland
Abstract. Starting from a biochemical signalling pathway model ex-
pressed in a process algebra enriched with quantitative information we
automatically derive both continuous-space and discrete-state represen-
tations suitable for numerical evaluation. We compare results obtained
using implicit numerical differentiation formulae to those obtained using
approximate stochastic simulation thereby exposing a flaw in the use of
the differentiation procedure producing misleading results.
1 Introduction
The malfunction of cellular signalling processes has significant detrimental ef-
fects, leading to uncontrolled cell proliferation, as in cancer; or leading to other
cells in the body being attacked, as in auto-immune diseases. The dynamics of
cell signalling mechanisms are profoundly complex and at present are not fully
understood. Computational modelling of cell signal transduction is an important
intellectual tool in the scientific study of the biological processes which control
and regulate cellular function.
An example of an influential computational study of intracellular signal net-
works is [1]. The authors develop an ordinary differential equation (ODE) model
of epidermal growth factor (EGF) receptor signal pathways in order to give in-
sight into the activation of the MAP kinase cascade through the kinases Raf,
MEK and ERK-1/2. The ODE model is substantial, consisting of 94 state vari-
ables and 95 parameters. It is analysed using the numerical integration proce-
dures of the Matlab numerical computing platform and tested using sensitivity
analysis. The results increase our understanding of EGF receptor signal trans-
duction and suggest avenues for experimental work to test hypotheses generated
from the computational model. Published in 2002 the article is highly regarded
and has subsequently been cited by as many as 150 other research papers.
We have previously proposed a method of investigating cell signalling path-
ways using a process algebra enhanced with quantitative information, PEPA [2],
applied in [3] and [4]. Process algebras are well-known in theoretical computer
science but are still unfamiliar to most computational biologists so we wished to
C. Priami (Ed.): CMSB 2006, LNBI 4210, pp. 63–77, 2006.
c
Springer-Verlag Berlin Heidelberg 2006
64 M. Calder et al.
PEPA
ODE
dy
dt
'
&
$
%
SSA
P(τ, μ)dτ
'
&
$
%
@
@
@R
Fig. 1. A high-level model in the PEPA process algebra can be used to generate either
a system of ODEs or a stochastic simulation
help to establish their relevance by reproducing the results of [1], starting from
the published paper together with its supplementary material and the Matlab
ODE model made available by the authors.
We were able to reproduce the results from [1] starting from our model in
the PEPA process algebra but because we were starting from the vantage point
of modelling in process algebra we could apply other analysis procedures, un-
available to the authors of [1] (Figure 1). To our surprise when modelling in
process algebra we discovered that the computational simulation conducted by
ODEs in [1] contains a systematic flaw in the analysis process which affects many
of the results, some significantly. To the best of our knowledge these errors are
presently unknown: at the very least they were unknown to us. Using the insights
obtained from our analysis procedures we were able to return to the differen-
tial equation model, diagnose and correct the flaws in the analysis, and show
agreement between the results obtained using continuous-space analysis and the
results obtained using a discrete-state stochastic analysis.
Computational methods are well-understood to be complex and delicate so the
relevance of this finding is not that there is an error in one particularly rich and
valuable numerical study, or that modelling with ODEs is an unsatisfactory pro-
cedure, but rather that modelling in high-level languages (such as process algebras
or Petri nets) may give a methodological advantage which allows an entire class of
hard-to-detect errors and corner cases to be discovered and diagnosed before the
results are published and promulgated to the wider scientific community.
As original contributions the present paper contains the analysis of the process
used to detect the error in the earlier modelling study [1], a description of the new
software tool used for integrated continuous-space and discrete-state stochastic
analysis of PEPA process algebra models, and an overview of an extensive pro-
cess algebra modelling study comprising 188 process definitions describing the
dynamics of 95 of the reaction channels in the signalling cascade of the EGF
receptor-induced MAP kinase pathway.
Structure of this paper: In Section 2 we present background material on our
previous work. We follow this in Section 3 with a discussion of related work. In
Section 4 we present an introduction to quantitative process algebras, considering
Stronger Computational Modelling of Signalling Pathways 65
the expressive capabilities of these languages. In Section 5 we explain how these
languages are used in modelling. Section 6 presents a comparison of our analysis
results and the results of other authors. In Section 7 we discuss the software tool
used to perform the analysis. Finally, we present conclusions in Section 8.
2 Background
In an earlier study we made two distinct computational models of the Ras/Raf-
1/MEK/ERK signalling pathway, both expressed in the PEPA process algebra.
Our models were based on the deterministic model presented directly as a system
of coupled ordinary differential equations in [5].
Our process algebra models adhere to two distinct modelling styles—the
reagent-centric and pathway models from [3]. We interpreted these under the
continuous-time Markov chain semantics for the PEPA language, and thus these
gave rise to stochastic models of the pathway. We used well-known procedures
of numerical linear algebra to conduct a quantitative stochastic evaluation of
the pathway. We used the process algebraic reasoning apparatus of the PEPA
language to establish that these two models were strongly equivalent, meaning
that a timing-aware observer could not distinguish between them. In the exten-
sion of this work in [6] we presented automatic procedures for converting in both
directions between the reagent-centric and pathway views.
We revisited the reagent-centric model in [4], mapping it to a system of ODEs.
The model considered in [4] adds additional species to the model presented
in [5] in order to concentrate on a detail of the pathway not considered in [5].
We applied the mapping procedure from [4] to a reduced version of the model
without these additional species and were able to show that the model gave rise
to exactly the same system of ODEs as studied previously in [5] establishing
a precise formal equivalence between the process algebra model and the ODE
model.
The deterministic and stochastic approaches to computational modelling in
systems biology are often presented as alternatives; one should choose one ap-
proach or the other. Some authors have suggested that stochastic approaches are
technically superior because they can expose small-scale effects which are caused
by some molecular species being present in the reaction volume in very low copy
numbers. We are instead in agreement with the authors of [7], who argue that
the principal challenge is choosing the appropriate framework for the modelling
study at hand. For some problems the influence of effects such as intra-cellular
noise or circumstances such as low copy numbers is sufficiently great that a
thorough stochastic treatment is essential. In other modelling problems no such
influences are manifest and a deterministic treatment based on reaction rate
equations is the correct approach.
The divergence between the stochastic behaviour exposed at low copy num-
bers of reactants and the deterministic approach based on reaction rate equations
is due to the reliance of the ODE-based analysis on the assumption of continuity
and the use of the law of mass action, essentially an empirical law derived from
66 M. Calder et al.
in vitro experimentation. Gillespie’s Stochastic Simulation Algorithm (SSA) [8]
makes no use of such an empirical law, and is instead grounded in the theory of
statistical thermodynamics. In consequence it is an exact procedure for numer-
ically simulating the dynamic evolution of a chemically reacting system, even
at low copy numbers. However, the SSA method converges, as the number of
reactants increases, to the solution computed by the ODEs so that the methods
are in agreement in the limit [9].
Gillespie’s exact algorithm models systems in which there are M possible
reactions represented by the indexed family R
µ
(1 ≤ μ ≤ M ). It builds on a
reaction probability density function P (τ,μ | X) such that P (τ,μ | X)dτ is the
probability that given the state X at time t,thenext reaction in the volume will
occur in the infinitesimal time interval (t + τ, t+ τ +dτ) and be an R
µ
reaction.
Starting from an initial state, SSA randomly picks the time and type of the next
reaction to occur, updates the global state to record the fact that this reaction
has happened, and then repeats.
In practice, Gillespie’s SSA is effective only for non-stiff systems on short
time scales. An approximate acceleration procedure called “τ-leaping” was later
developed by Gillespie and Petzold [10]. The “implicit τ-leaping” method [11]
was developed to attack the orthogonal problem of stiffness, common in multi-
scale modelling, where different time-scales are appropriate for reactions. Recent
advances in the field include the development of slow-scale SSA which produces
a dramatic speed-up relative to SSA by prioritising rare events [12].
A recent survey paper on stochastic simulation is [13]. A comparison paper on
stochastic simulation methods and their relation to differential-equation based
analysis of reaction kinetics is [9].
3 Related Work
We are not the first authors to investigate the model from [1] using stochastic
simulation methods. An earlier comparison using the binomial τ-leap method
appeared in [14]. However, the authors of [14] compare the solutions computed
by their binomial τ-leap method with the solutions computed by Gillespie’s
stochastic simulation algorithm and did not compare with the results from [1].
For this reason the authors of [14] did not find the error which we uncovered
by comparing the results computed by stochastic simulation with the results
computed by the authors of [1] using ordinary differential equations.
In [15] the authors use the PRISM probabilistic model checker [16] to check
logical formulae of Continuous Stochastic Logic (CSL) [17] against models of sig-
nalling pathways expressed as state-machines in the PRISM modelling language,
comparing the result against an ODE model coded in the Matlab numerical
platform.
A recent technical note [18] uses modelling in a stochastic process calculus and
stochastic simulation to investigate the MAPK cascade previously studied in [19]
using ordinary differential equations. [18] uses synthetic values for rate constants
(all are set to 1.0) so comparison with the results of [19] is not meaningful.
Stronger Computational Modelling of Signalling Pathways 67
4ProcessAlgebras
Process algebras are concise formally-defined modelling languages for the precise
description of concurrent, communicating systems. Our belief is that they are
well-suited to modelling cell signalling pathways and our interest here is exclu-
sively in process algebras which are decorated with quantitative information [20].
The PEPA process algebra [2] which we use benefits from formal semantic de-
scriptions of different characters which are appropriate for different uses. The
structured operational semantics presented in [2] maps the PEPA language to
a Continuous-Time Markov Chain (CTMC) representation. A continuous-space
semantics maps PEPA models to a system of ordinary differential equations
(ODEs) [21], admitting different solution procedures.
4.1 Expressiveness
Because we are modelling in a high-level language it is possible to apply these
very different numerical evaluation procedures to compute different kinds of
quantitative information from the same model. This is a freedom which we would
not have if we had coded a Markov chain or a differential equation-based rep-
resentation of the model directly in a numerical computing platform such as
Matlab. One freedom which the use of a high-level language gives the modeller
is the possibility to use either discrete-state or continuous-space analysis pro-
cedures. Another is the option of applying both types of analysis to the same
model, and that is the approach which we have used here.
One strength of the PEPA process algebra as an expressive and practical
modelling language is its support for multi-way co-operation;wehavemade
use of this expressive power in all of our modelling studies in systems biology.
Genuinely tri-molecular collisions occur only exceptionally rarely in dilute fluids
so these do not normally arise in our modelling for this reason. Rather a collision
between, say, an enzyme and a substrate to produce a compound, is expressed
in PEPA as a three-way co-operation between the input enzyme and substrate
(whose molecular concentrations are reduced) and the output compound (whose
molecular concentration is increased). Similarly a reaction channel with two
input species and two output species is represented as a four-way co-operation
in PEPA. Some reaction channels may have more inputs or more outputs and
so having this expressive power available in our chosen process algebra seems
well-suited to the type of modelling which is undertaken in the area.
4.2 Combinators of the Language
We give only a brief introduction to the PEPA language here. The reader is
referred to [2] for the definitive description.
PEPA provides a set of combinators which allow expressions to be built which
define the behaviour of components via the activities that they engage in. These
combinators are presented below.
![Fig. 5. Graphs of differential equation and stochastic simulation results compared. The solid red line is the solution of the original model from [1], which shows marked differences in some graphs from the solution of the PEPA-derived τ -leap simulation and (a dashed green line) and the solution of the ODE model using smaller time steps (a dotted blue line). The solution of the PEPA-derived τ -leap simulation and the solution of the ODE model using smaller time steps are virtually indistinguishable in the graphs.](/figures/fig-5-graphs-of-differential-equation-and-stochastic-2vp4k199.webp)
