Stochastic Simulation of Process Calculi for Biology

Abstract

Biological systems typically involve large numbers of components with complex, highly parallel interactions and intrinsic stochasticity. To model this complexity, numerous programming languages based on process calculi have been developed, many of which are expressive enough to generate unbounded numbers of molecular species and reactions. As a result of this expressiveness, such calculi cannot rely on standard reaction-based simulation methods, which require fixed numbers of species and reactions. Rather than implementing custom stochastic simulation algorithms for each process calculus, we propose to use a generic abstract machine that can be instantiated to a range of process calculi and a range of reaction-based simulation algorithms. The abstract machine functions as a just-in-time compiler, which dynamically updates the set of possible reactions and chooses the next reaction in an iterative cycle. In this short paper we give a brief summary of the generic abstract machine, and show how it can be instantiated with the stochastic simulation algorithm known as Gillespie’s Direct Method. We also discuss the wider implications of such an abstract machine, and outline how it can be used to simulate multiple calculi simultaneously within a common framework.

1 Introduction

Biological systems typically involve large numbers of components with complex, highly parallel interactions and intrinsic stochasticity. To model this complexity, numerous programming languages based on process calculi have been developed, many of which are expressive enough to generate unbounded numbers of molecular species and reactions. Examples include variants of the stochastic pi-calculus [14, 15, 12], BlenX [4], the kappa calculus [3], and variants of the bioambient calculus [16, 11]. As a result of this expressiveness, such calculi cannot rely on standard reaction-based simulation methods such as [8, 6], which require fixed numbers of species and reactions. Instead, a custom simulation algorithm is typically developed for each calculus. The choice of algorithm depends on the nature of the underlying biological system, such as whether exact simulation is required [7, 6], whether certain reactions operate at different timescales [8, 17], or whether non-Markovian reaction rates are needed [2, 9].

Rather than implementing custom stochastic simulation algorithms for each process calculus, we propose to use a generic abstract machine that can be instantiated to a range of process calculi and a range of reaction-based simulation algorithms. The abstract machine functions as a just-in-time compiler, which dynamically updates the set of possible reactions and chooses the next reaction in an iterative cycle. The abstract machine is instantiated to a particular calculus by defining two functions: one for transforming a process of the calculus to a set of species, and another for computing the set of possible reactions between species. The abstract machine is instantiated to a particular simulation algorithm by definition three functions: one for computing the next reaction, one for computing the reaction activity from an initial set of reactions and species populations, and a third for updating the reaction activity as the species populations change over time. Having a clear separation between the simulation algorithm and the language specification allows us not only to easily instantiate the machine to different process calculi, but also to add new simulation algorithms that can be shared between calculi. Furthermore, the approach could be used to dynamically integrate the simulation of multiple process calculi simultaneously, acting as a common language runtime for the simulation of process calculi for biology.

In this short paper we give a brief summary of the generic abstract machine of [9], and show how it can be instantiated with the stochastic simulation algorithm of [7]. We also discuss the wider implications of such an abstract machine, and outline how it can be used to simulate multiple calculi simultaneously within a common framework.

2 Summary of the Abstract Machine

Table 1: Syntax of the generic abstract machine, where a term $T$ consists of the current time $t$, a species map $S$ and a reaction map $R$. We let $\widetilde{I}$ denote a multiset of species $\{I_{1},..,I_{N}\}$ and $\widetilde{O}$ denote a set of reactions.

T	syntax	description
$T$	$(t,S,R)$	Time $t$, species map $S$, reaction map $R$
$S$	$\{I_{1}\mapsto i_{1},..,I_{N}\mapsto i_{N}\}$	Map from a species $I$ to its population $i$,
$R$	$\{O_{1}\mapsto A_{1},..,O_{N}\mapsto A_{N}\}$	Map from a reaction $O$ to its activity $A$
$O$	$(\widetilde{I},r,\widetilde{I^{\prime}})$	Reaction with reactants $\widetilde{I}$, products $\widetilde{I^{\prime}}$ and rate $r$.

Table 2: Parameterised definition of the generic abstract machine. If $\widetilde{I}$ is a multiset $\{I_{1},..,I_{N}\}$ we write $\widetilde{I}\oplus T$ for $I_{1}\oplus..\oplus I_{N}\oplus T$, and $T\ominus\widetilde{I}$ for $T\ominus I_{1}\ominus..\ominus I_{N}$ (the order is unimportant). We write ${\it dom}(S)$ for the domain of $S$. We also write $S(I)$ for the value associated with $I$ in $S$, and $S\{I\mapsto v\}$ for $S$ updated so that $v$ is associated with $I$.

function		definition
$P\oplus T$	$\triangleq$	${\it species}(P)\oplus T$
$I\oplus(t,S,R)$	$\triangleq$	$(t,S^{\prime},R\cup R^{\prime})$ if $\widetilde{I^{\prime}}={\it dom}(S)$; $I\notin\widetilde{I^{\prime}}$; $S^{\prime}=S\{I\mapsto 1\}$;
		$\widetilde{O}={\it reactions}(I,\widetilde{I^{\prime}})$; $R^{\prime}={\it init}(\widetilde{O},(t,S^{\prime},R))$
$I\oplus(t,S,R)$	$\triangleq$	$(t,S^{\prime},R\cup R^{\prime})$ if $S(I)=i$; $S^{\prime}=S\{I\mapsto i+1\}$; $R^{\prime}={\it updates}(I,(t,S^{\prime},R))$
$(t,S,R)\ominus I$	$\triangleq$	$(t,S^{\prime},R\cup R^{\prime})$ if $S(I)=i$; $S^{\prime}=S\{I\mapsto i-1\}$; $R^{\prime}={\it updates}(I,(t,S^{\prime},R))$

The syntax of the generic abstract machine is summarised in Table 2, and is based on the definitions of [9]. A machine term $T$ is a triple $(t,S,R)$, where $t$ is the current time, $S$ is a map from a species $I$ to its population $i$, and $R$ is a map from a reaction $O$ to its activity $A$, which is used to compute the next reaction. Each reaction is represented by a tuple $(\widetilde{I},r,\widetilde{I^{\prime}})$, where $\widetilde{I}$ denotes the multiset of reactant species, $\widetilde{I^{\prime}}$ denotes the multiset of product species and $r$ denotes the reaction rate. The structure of a term of the abstract machine can be summarised as follows.

Machine term $T$
Time $t$	Species map $S$		Reaction map $R$
	Species	Population	Reaction	Activity
	$I_{1}$	$i_{1}$	$\widetilde{I}_{1}\stackrel{{\scriptstyle r_{1}}}{{\longrightarrow}}\widetilde{I}_{1}^{\prime}$	$A_{1}$
	$\ldots$	$\ldots$	$\ldots$	$\ldots$
	$I_{N}$	$i_{N}$	$\widetilde{I}_{M}\stackrel{{\scriptstyle r_{M}}}{{\longrightarrow}}\widetilde{I}_{M}^{\prime}$	$A_{M}$

To instantiate the abstract machine with a given process calculus, it is sufficient to define a function ${\it species}(P)$ for transforming a process $P$ to a multiset of species, together with a function ${\it reactions}(I,\widetilde{I^{\prime}})$ for computing the set of reactions between a new species $I$ and an existing set of species $\widetilde{I^{\prime}}$. The syntax of species $I$ is specific to the choice of process calculus. The ${\it species}$ function is used to initialise the abstract machine at the beginning of a simulation, while the ${\it reactions}$ function is used to update the set of possible reactions dynamically. This allows systems with potentially unbounded numbers of species and reactions to be simulated.

To instantiate the abstract machine with a given simulation algorithm, it is sufficient to define a function ${\it next}(T)$ for choosing the next reaction from a term $T$, a function ${\it init}(\widetilde{O},T)$ for initialising a term with a set of reactions $\widetilde{O}$, and a function ${\it updates}(I,T)$ for updating the reactions in a term affected by a given species $I$. The abstract machine is then executed by repeated application of the following rule.

$$\dfrac{(\widetilde{I},r,\widetilde{I^{\prime}}),t^{\prime}={\it next}(t,S,R)}{t,S,R\stackrel{{\scriptstyle(\widetilde{I},r,\widetilde{I^{\prime}})}}{{\longrightarrow}}\widetilde{I^{\prime}}\oplus((t^{\prime},S,R)\ominus\widetilde{I})}$$

Each time the next reaction is selected, it is executed by removing the reactants $\widetilde{I}$ from the machine term, adding the products $\widetilde{I^{\prime}}$ and updating the current time of the machine. Corresponding definitions for adding and removing species are summarised in Table 2. A process $P$ is added to a machine term $T$ by computing the multiset of species $\{I_{1},\dots,I_{N}\}$ which correspond to $P$ and then adding each of these species to the term. If a new species $I$ is already present in the term then its population is incremented in $S$ and the activity of the affected reactions is updated. If the species is not already present in the term, its population is set to $1$ in $S$ and new reactions for the species are computed, together with their activity. The operation $T\ominus\widetilde{I}$ removes the species $\widetilde{I}$ from the machine term $T$, by decrementing the corresponding species populations and by updating the affected reactions.

3 Instantiating the Abstract Machine

Table 3: Instantiation of the generic abstract machine with the stochastic simulation algorithm of [7]. We write $\{E_{i}\mid C_{1};..;C_{N}\}$ to denote the set of elements $E_{i}$ that satisfy conditions $C_{1};..;C_{N}$. We let $n_{1}$ and $n_{2}$ denote two random numbers from the standard uniform distribution, $U(0,1)$. The function ${\it propensity}(O,S)$ is defined in the main text.

function		definition
$\mathit{next}(t,S,R)$	$\!\triangleq\!$	$O_{\mu},t\!+t^{\prime}$ if $a_{0}=\sum_{O_{i}\in{\it dom}(R)}R(O_{i})$; $t^{\prime}=(\frac{1}{a_{0}})ln(\frac{1}{n_{1}})$; $\sum_{i=1}^{\mu-1}a_{i}<n_{2}a_{0}\leq\sum_{i=1}^{\mu}a_{i}$
$\mathit{init}(\widetilde{O},(t,S,R))$	$\!\triangleq\!$	$\{O_{i}\mapsto{\it propensity}(O_{i},S)\>\mid\>O_{i}\in\widetilde{O}\}$
$\mathit{updates}(I,(t,S,R))$	$\!\triangleq\!$	$\{O_{i}\mapsto{\it propensity}(O_{i},S)\>\mid\>O_{i}\in{\it dom}(R);O_{i}=(\tilde{J},r,\tilde{J^{\prime}});I\in\tilde{J}\}$

An instantiation of the abstract machine with the stochastic simulation algorithm of [7] is outlined in Table 3. Each reaction $(\widetilde{I},r,\widetilde{I^{\prime}})$ is mapped to its propensity $a_{i}$, which is computed by multiplying the rate of the reaction by the number of distinct combinations of the reactants $\widetilde{I}$. The function ${\it propensity}(O,S)$ computes the propensity of the reaction $O$ given the species map $S$ and is defined as follows, assuming that reactions are either unary or binary and that $I_{1}$ and $I_{2}$ are distinct species.

$${\it propensity}((\{I_{1}\},r,\widetilde{I^{\prime}}),S) \triangleq r\times S(I_{1})$$

$${\it propensity}((\{I_{1},I_{1}\},r,\widetilde{I^{\prime}}),S) \triangleq r\times S(I_{1})\times(S(I_{1})-1)/2$$

$${\it propensity}((\{I_{1},I_{2}\},r,\widetilde{I^{\prime}}),S) \triangleq r\times S(I_{1})\times S(I_{2})$$

The function ${\it init}(\widetilde{O},T)$ computes the initial propensity for each reaction in $\widetilde{O}$, using the initial species populations in $T$, while the function ${\it updates}(I,T)$ updates the propensities of all the reactions in $T$ for which $I$ is a reactant. Finally, the function ${\it next}(T)$ chooses a reaction from $T$ with probability proportional to the reaction propensity, and computes the corresponding duration of the reaction according to [7].

We have also instantiated the abstract machine to the Next Reaction Method of [6] and to the Non-Markovian Next Reaction Method of [9], by defining corresponding ${\it init}$, ${\it next}$ and ${\it updates}$ functions, as described in [9]. We have used the abstract machine to implement the DNA Strand Displacement (DSD) calculus for modelling DNA circuits [13], the Genetic Engineering of Cells (GEC) calculus for modelling of genetic devices [10], and the Stochastic Pi Machine (SPiM) calculus for general modelling of biological systems [18], by defining appropriate $species$ and ${\it reactions}$ functions for each calculus. Simulators for these three calculi are available online at http://research.microsoft.com/dna, http://research.microsoft.com/gec and http://research.microsoft.com/spim, respectively. Technical details of the instantiation of the generic abstract machine with the stochastic pi-calculus and the bioambient calculus are outlined in [9]. We are currently developing an instantiation of the generic abstract machine to the kappa calculus of [3]. Although the idea of integrating different modelling and simulation methods within a common framework is not a new one [5], our approach is the first attempt to formally define a generic framework for simulating a broad range of process calculi with an arbitrary reaction-based simulation algorithm.

4 A Common Simulation Framework

The generic abstract machine can be used to simulate multiple calculi simultaneously by assuming a separate species type $I_{L}$ for each calculus $L$, together with an initial set of cross-calculus reactions $\widetilde{O_{0}}$. An example of a cross-calculus reaction is $I_{DSD}+I_{SPiM}\stackrel{{\scriptstyle r}}{{\longrightarrow}}I_{SPiM}+I^{\prime}_{SPiM}$, which takes a species of the DSD language, such as a known DNA vaccine assembled via strand displacement, together with a species of the SPiM language, such as a polymerase, and produces a corresponding protein species in SPiM together with the original polymerase. The reaction therefore enables the output of a strand displacement model in DSD to interface with a cellular model in SPiM. For each dynamically created species $I_{L}$ the function ${\it reactions}(I_{L},\widetilde{I^{\prime}})$ calls the appropriate calculus-specific function ${\it reactions_{L}}(I_{L},\widetilde{I_{L}^{\prime}}),$ where $\widetilde{I_{L}^{\prime}}$ denotes the subset of species in $\widetilde{I^{\prime}}$ that are of type $L$. This approach allows multiple calculi to interact with each other within the same simulation environment, via a fixed set of interface reactions. Further work is needed to formalise the multi-language execution paradigm in more detail.

The generic abstract machine can therefore be used to simulate a range of existing process calculi within a common framework. By decoupling the choice of calculus from the choice of simulation algorithm, multiple calculi can re-use the same algorithm via a common interface, without the need to implement custom simulation algorithms for each calculus. In future, this could allow models to be constructed from components written in different domain-specific languages, each designed to allow a natural, concise encoding of that component. The components could then interact dynamically via a common language runtime, allowing integrated simulation of heterogeneous biological systems.