Implementing the Stochastics Brane Calculus
in a Generic Stochastic Abstract Machine
Marino Miculan

Ilaria Sambarino

Department of Mathematics and Computer Science, University of Udine, Italy
marino.miculan@uniud.it

ilaria.sambarino@gmail.com

In this paper, we deal with the problem of implementing an abstract machine for a stochastic version of the Brane Calculus. Instead of defining an ad hoc abstract machine, we consider the generic
stochastic abstract machine introduced by Lakin, Paulevé and Phillips. The nested structure of membranes is flattened into a set of species where the hierarchical structure is represented by means of
names. In order to reduce the overhead introduced by this encoding, we modify the machine by
adding a copy-on-write optimization strategy. We prove that this implementation is adequate with respect to the stochastic structural operational semantics recently given for the Brane Calculus. These
techniques can be ported also to other stochastic calculi dealing with nested structures.

1

Introduction

A fundamental issue in Systems Biology is modelling the membrane interaction machinery. Several
models have been proposed in the literature [11, 14, 3]; among them, the Brane Calculus [4] has been
arisen as a good model focusing on abstract membrane interactions, still being sound with respect to
biological constraints (e.g. bitonality). In this calculus, a process represents a system of nested compartments, where active components are on membranes, not inside them. This reflects the biological
evidence that functional molecules (proteins) are embedded in membranes, with consistent orientation.
In the original definition of the Brane Calculus [4] (which we will recall in Section 2) membranes interact according to three basic reaction rules corresponding to phagocytosis, endo/exocytosis, and pinocytosis. However, this semantics does not take into account quantitative aspects, like stochastic distributions, which are important for, e.g., implementing stochastic simulations.
A stochastic semantics for the Brane Calculus has been provided in [2], following an approach pioa,r
neered in [5] (but see also [8, 12] for Markov processes). Instead of giving a stochastic version P −→ Q
of the reaction relation, in this semantics each process is given a measure of the stochastic distribution
of the possible outcomes. More precisely, we define a relation P → µ associating to a process P an
action-indexed family of measures µ: for an action a, the measure µa specifies for each measurable set
S of processes, the rate µa (S) ∈ R+ of a-transitions from P to (elements of) S. An advantage of this
approach is that we can apply results from measure theory for solving otherwise difficult issues, like
instance-counting problems; moreover, process measures are defined compositionally, and in fact the
relation P → µ can be characterized by means of a set of rules in a GSOS-like format. We will recall this
stochastic semantics and its main properties in Section 3.
In this paper, we use this new semantics for defining a stochastic abstract machine for the Brane
Calculus, so that it can be effectively used for in silico simulations of membrane systems. Defining an
ad hoc abstract machine for the Brane Calculus would be a complex task; instead, we take advantage of
the generic abstract machine for stochastic process calculi (GSAM for short) introduced in [13, 10] as
a general tool for simulating a broad range of calculi. This machine can be instantiated to a particular
Gabriel Ciobanu (Ed.): Membrane Computing and
Biologically Inspired Process Calculi (MeCBIC 2012)
EPTCS 100, 2012, pp. 82–100, doi:10.4204/EPTCS.100.6

83

M. Miculan and I. Sambarino

calculus by defining a function for transforming a process of the calculus to a set of species, and another
for computing the set of possible reactions between species.
An important aspect is that this abstract machine does not have a native notion of compartment,
which is central in the Brane Calculus (as in any other model of membranes). To overcome this problem,
we adopt a “flat” representation of membrane systems, used also in [10], where the hierarchical structure
is represented by means of names: each name represents a compartment, and each species is labelled
with the name of the compartment where it is located, and the name of its inner compartment (if any).
So names and species are the nodes and the arcs of the tree, respectively. This technique can be used for
representing any system with a tree-like structure of compartments.
However, this approach does not scale well, as the population of species may grow enormously: for
instance, a population of n identical cells would lead to n species, all differing only for the name of its
inner compartment, instead of a single specie with multiplicity n. For circumventing this problem, in
Section 4 we introduce a variant of the GSAM with a copy-on-write optimization strategy—hence called
COWGSAM. The idea is to keep a single copy of each species, with its multiplicity; when a reaction has
to be applied, fresh copies of the compartments involved are generated on-the-fly, and reactions and rates
are updated accordingly. In this way, the hierarchical structure is unfolded only if and when needed.
In Section 5 we show how the Brane Calculus can be represented in the COWGSAM, and we will
prove that the abstract machine obtained in this way is adequate with respect to the stochastic semantics
of the Brane Calculus; in this proof, we take advantage of the compositional definition of this semantics.
Conclusions and final remarks are in Section 6.

2

Brane Calculus

In this section we recall Cardelli’s Brane Calculus [4] focusing on its basic version (without communication primitives, complexes and replication).
First, let us fix the notation we will use hereafter. Let S be a set of sorts (or “types”), ranged over
by s,t, and T a set of S-sorted terms; for t ∈ S, Tt ⊆ T denotes the set of terms of sort t. For A a set of
symbols, A∗ denotes the set of finite words (or lists) over A, and ha1 , . . . , an i denotes a word in A∗ . For a
word ht1 , . . . ,tn i in S∗ , we define Tht1 ,...,tn i , Tt1 × · · · × Ttn .
Syntax

The sorts and the set B of terms of Brane Calculus are the following:

Sorts :: S
Membranes :: Bmem
Systems :: Bsys

t ::= sys | mem

σ , τ ::= 0 | σ |τ | Jn .σ | JIn (τ).σ | Kn .σ | KIn .σ | Gn (τ).σ
P, Q ::= k | P m Q | σ hPi

The subscripted names n are taken from a countable set Λ. By convention we shall use M, N, . . . to
denote generic Brane Calculus terms in B.
A membrane can be either the empty membrane 0, or the parallel composition of two membranes
σ |τ, or the action-prefixed membrane ε.σ . Actions are: phagocytosis J, exocytosis K, and pinocytosis
G. Each action but pinocytosis comes with a matching co-action, indicated by the superscript ⊥ .
A system can be either the empty system k, or the parallel composition P m Q, or the system nested
within a membrane σ hPi. Notice that, differently from [4], pino actions are indexed by names in Λ.
In [4], names are meant only to pair up an action with its corresponding co-action, hence a pino action
does not need to be indexed by any name. Actually, names can be thought of as an abstract representation

84

Implementing Stochastics Brane Calculus in COWGSAM

JIn (ρ).τ|τ0 hQi m Jn .σ |σ0 hPi } τ|τ0 hρhσ |σ0 hPii m Qi
KIn .τ|τ0 hKn .σ |σ0 hPi m Qi } σ |σ0 |τ|τ0 hQi m P
G(ρ).σ |σ0 hPi } σ |σ0 hρhki m Pi
P}Q
PmR}QmR

(red-comp)

(red-exo)

P}Q
σ hPi } σ hQi

(red-pino)

P ≡ P0

(red-phago)

P0 } Q0
P}Q

Q0 ≡ Q

(red-loc)

(red-equiv)

Table 1: Reduction semantics for the Brane Calculus.

of particular protein conformational shapes; hence, each name can correspond to a different biological
behaviour. Therefore, if we want to observe also kinetic properties of processes, it is important to keep
track of names in pino actions.
Terms can be rearranged according to a structural congruence relation; the intended meaning is
that two congruent terms actually denote the same system. Structural congruence ≡ is the smallest
equivalence relation over B which satisfies the axioms and rules listed below.
PmQ ≡ QmP

σ |τ ≡ τ|σ

0hki ≡ k

P m (Q m R) ≡ (P m Q) m R
σ |(τ|ρ) ≡ (σ |τ)|ρ

P≡Q
PmR ≡ QmR

α ∈ {Jn , Kn , KIn }n∈Λ
α.σ ≡ α.τ

σ ≡τ

σ ≡τ
σ |ρ ≡ τ|ρ

Pmk ≡ P

σ |0 ≡ σ

P≡Q σ ≡τ
σ hPi ≡ τhQi

β ∈ {JIn , Gn }n∈Λ ρ ≡ ν σ ≡ τ
β (ρ).σ ≡ β (ν).τ

Differently from [4], we allow to rearrange also the sub-membranes contained in co-phago and pino
actions (by means of the last inference rule above).
Reduction Semantics The dynamic behaviour of Brane Calculus is specified by means of a reduction
semantics, defined over a reduction relation (“reaction”) } ⊆ Bsys × Bsys , whose rules are listed in
Table 1. Notice that the presence of (red-phago/exo/pino) and (red-equiv) makes this not a structural
presentation, since these rules are not primitive recursive in the syntax (i.e., structural recursive) as
required by the SOS format.

3

Stochastic Structural Operational Semantics for the Brane Calculus

In this section we recall the stochastic structural operational semantics for the Brane Calculus, as defined
a,r
in [2]. Following [5], we replace the classic “pointwise” rules of the form P −→ P0 with rules of the form
P → µ, where µ is an indexed class of measures on the measurable space of processes. We assume the
reader to be familiar with basic notions from measure theory; for a brief summary, see Appendix A.

85

M. Miculan and I. Sambarino

0 →mem ωmem
β ∈ {JIn , Gn | n ∈ Λ}
β (τ).σ →mem [β ]τσ
k →sys ωsys

(void)

(zero)

(pref-arg)

α ∈ {Jn , Kn , KIn | n ∈ Λ}
α.σ →mem [α]σ

(pref)

σ →mem µ 0
τ →mem µ 00
σ |τ →mem µ 0 σ τ µ 00

σ →mem ν
P →sys µ
σ hPi →sys µ@σP ν

(loc)

(par)

P →sys µ 0
Q →sys µ 00
P m Q →sys µ 0 P ⊗Q µ 00

(comp)

Table 2: Stochastic structural operational semantics for Brane Calculus

The set of action labels for the Brane Calculus will be denoted by A and can be partitioned with
respect to the source sort (i.e., either systems or membranes), as follows:
Asys , {id : sys → hsysi} ∪ {phn : sys → hsys, sysi | n ∈ Λ} ∪
{ph⊥
n : sys → hmem, mem, sys, sysi | n ∈ Λ} ∪
{exn : sys → hmem, sys, sysi | n ∈ Λ}

Amem , {Jn , Kn , KIn : mem → hmemi | n ∈ Λ} ∪

{JIn , Gn : mem → hmem, memi | n ∈ Λ}

Let a range over A, and ar(a) denote its arity. To ease the reading in the following we will use the
N
notation ∆a (T, Σ) to denote the set of measures ∆(Tht1 ,...tn i , ni=1 Σti ), for ar(a) = t → ht1 , . . . ,tn i.
Let B/≡ be the set of ≡-equivalence classes on B. For M ∈ B, we denote by [M]≡ the ≡-equivalence
class of M (sometimes dropping the equivalence symbol when clear from the context).
Definition 3.1 (Measurable space of terms). The measurable space of terms (B, Π) is given by the
measurable space over B where Π is the σ -algebra generated by B/≡ .

Notice that B/≡ is a denumerable partition of B, hence it is a base (a generator such that all its
elements are disjoint) for Π. Any element of Π can be obtained by a countable union of elements of
S
the base, i.e., for all M ∈ Π there exist {Mi }i∈I , for some countable I, such that M = i∈I [Mi ]≡ . As a
consequence, in order to generate the whole Π we can simply compute all these unions, without the need
of any closure by complement.
N
A similar argument holds for the product space (Bht1 ,...,tn i , ni=1 Πti ), where ti ∈ {mem, sys} (1 ≤
Nn
i ≤ n); indeed i=1 Πti can be generated from the base Bht1 ,...,tn i /≡ht1 ,...,tn i , where ≡ht1 ,...,tn i ⊆ Bht1 ,...,tn i ×
Bht1 ,...,tn i is defined by
hM1 , . . . , Mn i ≡ht1 ,...,tn i hN1 , . . . , Nn i

iff

Mi ≡ Ni , for all 1 ≤ i ≤ n ,

which can be easily checked to be an equivalence relation. ≡ht1 ,...,tn i -equivalence classes are rectangles,
N
i.e. [hM1 , . . . , Mn i]≡ht1 ,...,tn i = [M1 ]≡ × · · · × [Mn ]≡ , therefore the product measure ni=1 Πti is well defined.
For sake of simplicity in the following we write [hM1 , . . . , Mn i]≡ in place of [hM1 , . . . , Mn i]≡ht1 ,...,tn i , and
Bht1 ,...,tn i /≡ in place of Bht1 ,...,tn i /≡ht1 ,...,tn i .
The operational semantics associates with each membrane a family of measures in ∆Amem (B, Π),
and with each system a family of measures in ∆Asys (B, Π). This can be represented by two relations

86

Implementing Stochastics Brane Calculus in COWGSAM

→mem : Tmem → ∆Amem (B, Π), →sys : Tsys → ∆Asys (B, Π), defined by the SOS rules listed in Table 2. (In
the following, for sake of readability, we will drop the indexes mem ,sys ). In these rules we use some
constants and operations over indexed families of measures, that we define next. For a set (of labels) A,
let us denote by ∆A (B, Π) the set ∏a∈A ∆a (B, Π) of A-indexed families of measures over (B, Π). Given a
family of measures µ ∈ ∆A (B, Π) and a ∈ A, the a-component of µ will be denoted as µa ∈ ∆a (B, Π).

Null: Let ωmem ∈ ∆Amem (B, Π) be the constantly zero measure, i.e., for all a∈Amem such that ar(a) =
N
t → ht1 , . . . ,tn i and M ∈ ni=1 Πti : (ωmem )a (M )=0.

Prefix: For arbitrary n ∈ Λ, α ∈ {J, K, KI }, and β ∈ {JI , G}, let the constants [αn ]σ , [βn ]τσ ∈ ∆Amem (B, Π)
be defined, for arbitrary X,Y ∈ Bmem /≡ , as follows:
(
ι(n) if n = m and σ ∈ X
([αn ]σ )αm (X) =
0
otherwise
([αn ]σ )βm (X ×Y ) = 0
([βn ]σ )αm (X) = 0
(
ι(n) if n = m and σ ∈ X, τ ∈ Y
([βn ]τσ )βm (X ×Y ) =
0
otherwise
Parallel: For µ, µ 0 ∈ ∆Amem (B, Π), let µ σ τ µ 0 ∈ ∆Amem (B, Π) be defined, for n ∈ Λ, α ∈ {J, K, KI },
S
β ∈ {JI , G}, and X,Y ∈ Bmem /≡ , as follows (where for all X, τ: X|τ , {[σ ]≡ | σ |τ ∈ X}):
(µ σ τ µ 0 )αn (X) = µαn (X|τ ) + µα0 n (X|σ )

(µ σ τ µ 0 )βn (X ×Y ) = (µ)βn (X|τ ×Y ) + (µ 0 )βn (X|σ ×Y )

Void: Let ωsys ∈ ∆Asys (B, Π) be defined by (ωsys )a (M ) = 0 for any a ∈ Asys , such that ar(a) = t →
N
ht1 , . . . ,tn i, and M ∈ ni=1 Πti .

Nesting: For ν ∈ ∆Amem (B, Π) and µ ∈ ∆Asys (B, Π), let µ@σP ν ∈ ∆Asys (B, Π) be defined, for X,Y ∈
Bmem /≡ and Z,W ∈ Bsys /≡ , as follows:
(
νJn ([σ ]≡ ) if σ hPi ∈ Z and k ∈ W
(µ@σP ν)phn (Z ×W ) =
0
otherwise
(
νJIn (X ×Y ) if P ∈ Z and k ∈ W
(µ@σP ν)ph⊥ (X ×Y × Z ×W ) =
n
0
otherwise
(
νKn (X) if P ∈ Z and k ∈ W
(µ@σP ν)exn (X × Z ×W ) =
0
otherwise

(µ@σP ν)id (X) = µid (Xσ hi ) +
n∈Λ

n∈Λ

∑

X 0 hX 00 h[k]≡ im[P]≡ i=X

νGn (X 0 × X 00 ) +

µexn (X 0 ×Y 0 ×Y 00 ) · νKIn (X 00 )
ι(n)
X 0 |X 00 hY 00 imY 0 =X

∑

87

M. Miculan and I. Sambarino

Composition: For µ, µ 0 ∈ ∆Asys (B, Π), let µ P ⊗Q µ 0 ∈ ∆Asys (B, Π) be defined, for X,Y ∈ Bmem /≡ and
S
Z,W ∈ Bsys /≡ , as follows (where for all W, Q, WmQ , {[P]≡ | P m Q ∈ W }):
0
(Z ×WmP )
(µ P ⊗Q µ 0 )phn (Z ×W ) = µphn (Z ×WmQ ) + µph
n
0
(µ P ⊗Q µ 0 )ph⊥ (X ×Y × Z ×W ) = µph⊥ (X ×Y × Z ×WmQ ) + µph
⊥ (X ×Y × Z ×WmP )
n

n

n

0
(µ P ⊗Q µ 0 )exn (X × Z ×W ) = µexn (X × Z ×WmQ ) + µex
(X × Z ×WmP )
n

(µ P ⊗Q µ 0 )id (X) = µid (XmQ ) + µid0 (XmP ) +
n∈Λ

∑

0 (X × X × Z × Z )
µphn (Y1 ×Y2 ) · µph
1
2
1
2
⊥
n

X1 hX2 hY1 imZ1 imY2 mZ2 =X
n∈Λ

∑

ι(n)

+

0 (Z × Z )
µph⊥ (X1 × X2 ×Y1 ×Y2 ) · µph
1
2
n
n

X1 hX2 hZ1 imY1 imZ2 mY2 =X

ι(n)

These operators have nice algebraic properties (e.g., µ 0 σ τ µ 00 = µ 00 τ σ µ 0 , (µ 0 σ τ µ 00 )σ |τ ρ µ 000 =
00
000
0
0
0
σ τ|ρ (µ τ ρ µ ), . . . ), and respect the structural congruence (e.g., if P ≡ P and Q ≡ Q then µ P ⊗Q
µ 00 = µ 0 P0 ⊗Q0 µ 00 ). We refer to [2] for further details about these properties, very useful in calculations.
The next lemmata state that the stochastic transition relation → (and hence the operational semantics)
is well-defined and consistent, that is, for each process we have exactly one family of measures of its
continuations, and this family respects structural congruence.
µ0

Lemma 3.2 (Uniqueness). For a ∈ A such that ar(a) = t → ht1 , . . . ,tn i, and M ∈ Bt , there exists a unique
µ ∈ ∆At (B, Π) such that M → µ.
Lemma 3.3. If M ≡ N and M → µ, then N → µ.

This operational semantics can be used to define the “ traditional” pointwise semantics:
a,r

4

M −→ hM1 , . . . , Mn i ⇐⇒ M → µ and µa ([hM1 , . . . , Mn i]≡ ) = r
and it is conservative with respect to the non-stochastic reduction semantics.
Proposition 3.4. For all systems P, Q ∈ Bmem , if P → µ and µid ([Q]) > 0 then P } Q.

4

The COW Generic Stochastic Abstract Machine

In this section we present a variant of the generic stochastic abstract machine (GSAM), oriented to
systems with nested compartments.
The GSAM has been introduced in [13, 10] for simulating a broad range of process calculi with an arbitrary reaction-based simulation algorithm. Although it does not have a native notion of “compartment”,
nested systems can represented by “flattening” all compartments and species into a single multiset, where
each species is tagged with names representing their position in the hierarchy, as shown in [10]. The idea
is to represent each compartment as a different species, keeping track of their position in the hierarchy
by means of (node) names. These names are ranged over by x, y, z, . . . , and are different from names in
actions. As an example, a system σ hτhii is represented as the multiset {σ hxy i 7→ 1, τhyz i 7→ 1}, which
means “there is one cell with membrane σ located in the compartment x and whose compartment is y,
and one cell with membrane τ positioned in y and whose compartment is z”. Reactions can happen only
if the names tagging the involved species match according to the required nesting structure.

88

Implementing Stochastics Brane Calculus in COWGSAM

Unfortunately, this approach does not scale well, as the population of species grows. Let us consider
a system composed of n copies of the same cell, e.g., n · (σ hi) (where n can be easily in the order of
103 –104 ). In the original GSAM idea, this should be represented in the machine as a single species with
multiplicity n, and each possible reaction is represented once but with propensity given by the law of mass
action taking into account the species’ multiplicity n. Instead, the “flat encoding” above yields n different
species σ hxy1 i, . . . , σ hxyn i, each with multiplicity 1; the set of reactions explodes correspondingly.
For circumventing this problem we introduce a variant of the GSAM with a copy-on-write strategy—
hence it is called COWGSAM. The idea is to keep a single copy of each species, with its multiplicity; the
same applies to reactions. When a reaction has to be applied, the compartments involved are “unfolded”,
i.e., fresh copies of the compartments are generated and the reaction set is modified accordingly; then,
the reaction can be applied. In this way, the hierarchical structure is unfolded only if and when needed.
In order to implement this idea, we have to modify the notion of machine term, reaction and reaction
rule. The COWGSAM (with the Next Reaction method) is shown in Figure 1.
First, for generating fresh names, we have to keep track of those already allocated. To this end we
introduce environments, which are finite sets of names. Then, the machine state is represented by a
machine term T , i.e. a quadruple E ` (t, S, R) where E is an environment; t is the current time; S is a
finite function mapping each species I to its population S(I) (if not null); and R maps each reaction O
to its activity A, which is used to compute the next reaction. (Notice that the syntax of species I is left
unspecified, as it depends on the specific process calculus one has to implement.) We say that a machine
term E ` (t, S, R) is well-formed if for all xi , x j ∈ E : xi = x j ⇒ i = j, and the free names in S, R appear in
E. In the following, we assume that machine terms are well formed.
r
Each reaction is a quadruple (S1 , r, f , S2 ), basically representing a reaction S1 →
− S2 , where
• S1 and S2 denote the reactant population and product population respectively;
• r denotes the rate (in s−1 ) of the reaction;

• f is a function mapping machine terms to machine terms; this functions implements any global
update of the machine term after the reaction (if needed).
a,O

The transitions of the abstract machine are represented by a relation T −−→ T 0 between machine
terms, indexed by the propensity a and reaction rules. This should be read as “T goes to T 0 with rate
a, using the rule O”. This relation is defined by (Reaction rule) in Figure 1, where the function next(T )
selects the next reaction, i.e. it returns a pair (O,t 0 ) where O = (S1 , r, f , S2 ) is the reaction to happen
first among all possible reactions in T , and t 0 is the new time of the system. Once the reaction has
been selected, we have first to create the separate (private) copies of the compartments involved, and to
update the reaction set accordingly. This is done by the functions cow( ) and dup( ), which implement
a deep copy-on-write: cow(E ` (t, S, R), S1 ) is a machine term E 0 ` (t, S0 , R0 ) representing the same state
as E ` (t, S, R), but in S0 the species indicated in S1 are unfolded; E 0 contains all names which have
been generated in the process, and R0 is the new set of reactions. (Actually, the freshly generated copies
represent the instances which are not involved by the reaction; this simplifies the reaction application.)
An example of the action of cow( ) is depicted in Figure 2.
At this point, the reaction is executed, by removing the reactants S1 from the machine term (via
the operation ), adding the products S2 (via the operation ⊕) and updating the current time of the
machine. The function f performs any global “clean-up” and restructuring to the machine term that
may be required by the reaction (e.g., garbage collection, elimination of names not used anymore,. . . ).
Moreover, since a reaction can rearrange the hierarchy structure, possibly creating new compartments
and deleting others, we have to add to the environment any fresh name introduced in the products.

89

M. Miculan and I. Sambarino

T ::=E ` (t, S, R)

(Machine term)

E ::=x1 , . . . , xn

(Environments)

S ::={I1 7→ i1 , . . . , IN 7→ iN }

(Populations)

R ::={O1 7→ A1 , . . . , ON 7→ AN }

(Reactions)

O ::=(S1 , r, f , S2 )

(Reaction)

((S1 , r, f , S2 ), a,t 0 ) = next(t, S, R) E 0 ` (t, S0 , R0 ) = cow(E ` (t, S, R), S1 )
a,(S1 ,r, f ,S2 )

E ` (t, S, R) −−−−−−→ norm(E 0 ∪ f n(S2 ) ` ( f (S2 ⊕ ((t 0 , S0 , R0 )

S1 ))))

(Reaction rule)

next(t, S, R) , (O, a,t 0 ) if R(O) = (t 0 , a) and t 0 = min{t | R(O) = (t, a)}
cow(E ` (t, S, R), 0)
/ , E ` (t, S, R)

cow(E ` (t, S, R), {ρhyz i 7→ j} ∪ S1 ) , cow(E 0 ∪ f n(S00 ) ` (t, S00 , R00 ∪ init(L0 ∪ L00 , (t, S0 , R)), S1 ) where
if S(σ hxy i) > 1 for some σ , x :
/E
let i = S(σ hxy i) and y0 ∈

(S0 , R0 ) = dup(E ` (t, S{σ hxy i 7→ 1, σ hxy0 i 7→ i − 1}, R), y, y0 )
L0 = reactions(σ hxy0 i 7→ i − 1, S)

otherwise let (S0 , R0 ) = (S, R), L0 = 0/ in
let n = S(ρhyz i), E 0 = E ∪ f n(S0 ), and z0 ∈
/ E0

(S00 , R00 ) = dup(E 0 ` (t, S0 {ρhyz i 7→ j, ρhyz0 i 7→ n− j}, R0 ), z, z0 )
L00 = reactions(ρhyz0 i 7→ n − j, S0 )

dup(E ` (t, S, R), y, y0 ) , (S0 ∪ S00 , R ∪ R00 ∪ init(L, (t, S0 , R)))
0

where S0 = S ∪ {ρhyw0 i 7→ i | S(σ hzy i) > 0, S(ρhyw i) = i}, w0 ∈
/ E,
0

L = reactions(ρhyw0 i 7→ i, S0 ),

(S00 , R00 ) = dup(E ∪ {w0 } ` (t, S, R), w, w0 )
Figure 1: The COW Generic Stochastic Abstract Machine, with the Next Reaction method.
Finally, the term can be “normalized” by collapsing equivalent copies of the same subtree into a
single copy (with the right multiplicity), by the function norm( ). In its simplest form, norm( ) can be
the identity, i.e., no normalization is performed at all. Although this is correct, it can lead to unnecessary
copies of the same subtrees. We can define norm( ) such that
norm({I1 7→ i1 , I2 7→ i2 } ∪ S, R) = norm({I1 7→ i1 + i2 } ∪ S, R[I2 /I1 ]) if I1 ≡ I2
where the equivalence between species can be implemented by comparing the subtrees starting from I1 ,
I2 , e.g., by calculating a suitable hash value. We leave this refinement as future development.
In order to implement the Next Reaction method, each reaction O is associated with a pair R(O) =
(a,t), where a is the reaction propensity and t is the time at which the reaction is supposed to occur. The
function delay(r, a) computes a time interval from a random variable with rate r and propensity a.

90

Implementing Stochastics Brane Calculus in COWGSAM

{I1 7→ i1 , . . . , IN 7→ iN } ⊕ (t, S, R) , I1 7→ i1 ⊕ . . . ⊕ IN 7→ iN ⊕ (t, S, R)

0
0
0
0
0

(t, S , R ∪ updates(I, (t, S , R))) if S(I) = i and S = S{I 7→ i + i}
I 7→ i ⊕ (t, S, R) , (t, S0 , R ∪ init(L, (t, S0 , R)))
if I ∈
/ dom(S), S0 = S{I 7→ i}


and L = reactions(I 7→ i, S)
(t, S, R)

{I1 7→ i1 , . . . , IN 7→ iN } , (t, S, R)
(t, S, R)

I1 7→ i1

...

IN 7→ iN

I 7→ i , (t, S0 , R ∪ updates(I, (t, S0 , R)))

if S(I) = i0 , i0 ≥ i and S0 = S{I 7→ i0 − i}

init(L, (t, S, R)) , {O 7→ (t 0 , a) | O ∈ L and a = propensity(O, S) and

O = (S1 , r, f , S2 ) and t 0 = t + delay(r, a)}

updates(I, (t, S, R)) , {O 7→ (t 0 , a0 ) | R(O) = (t 00 , a) and O = (S1 , r, f , S2 ) and S1 (I) > 0 and
if t 00 > t then a0 = propensity(O, S) and t 0 = t + (a/a0 )(t 00 − t)

if t 00 = t then a0 = propensity(O, S) and t 0 = t + delay(r, a)}
 ∗

 ∗

S (I1 )
S (In )
1 jn }
propensity((S1 , r, f , S2 ), S) , r ·
·...·
if S1 = {I1 7→ j1 , . . . , In 7→
j1
jn
(
S(σ hx i)
if x = root
S∗ (σ hxy i)(t0,, S0 , R0 ) ! yx(t1 , S1 , R∗ 1 ) !z · · · ! (tn , Sn , Rn ) ! . . .
S(σ hy i) · S (ρhx i) if x 6= root and S(ρhzx i) > 0
((S1 , r, f , S2 ), a,t
a,(S1 ,r, f ,S2 )

0)

= next(t, S, R)

Figure 1:EThe
COW Generic
Stochastic Abstract Machine (cont.).
` (t, S, R)
! E [ f n(S ) ` f (S ((t 0 , S, R) S ))
((S1 , r, f , S2

), a,t 0 )

E ` (t, S, R)
x
i·σ

= next(t, S, R)

a,(S1 ,r, f ,S2 )

!

E0

norm(E 0 [

2

`

2

(t, S0 , R0 )

1

= cow(E ` (t, S, R), S1 )

f n(S2 ) ` ( f (S2

((t 0 , S0 , R0 )

cow(E ` (t, S, R), S1 ) = E 0 ` (t, S0 , R0 )
cow(. . . , {rhyz i 7! j}})

S1 ))))

x

1·σ

(i-1)·σ

n·ρ y
j·ρ y (n-j)·ρ y' n·ρ
x
y
y
z x,y,w,z;w (P),
z' speciesx,y,w,z;z
z'' (Q)
speciesx;xz(s hthPi rhQii) = s hy i, thw i, rhz i, species
P speciesx (n · s hPi) = s hxy i, speciesy1 (P), . . . , s
Phxy i, species
P' yn (P) P''
n
1

= s hxy i 7! n, species
(P)
Figurespecies
2: Example
copy-on-write
of a ysubterm.
x (n · s hPi)
(S , r, f , S ) 2 R

1
2
This general structure can be instantiated with
a given
process calculus, just by providing the definition for the missing parts. Given a set Proc of processes,
we
have to define:
speciesx (n · s hm · thrhiii) = s hxy i 7! n, thyz i 7! m, rhzw i 7! 1
1. the syntax of the species I (which may be different from that of processes);

2. a function speciesE,x (P) mapping a process P ∈ Proc to a species set located in x;

speciesx (n·s hm·thrhiii) = s hx i 7! 1, thy1 i 7! 1, rhz11 i 7! 1, . . . , s hx i 7! 1, thy1 i 7! 1, rhznm i 7! 1

y1
z11
w11
yn
znm
wnm
3. a function reactions(I 7→ i, S) for
computing
the multiset
of reactions
between
a (new)
species I
with multiplicity i and a population of (existing) species S.

Jn .s hPi JIn (r).thQi } thQ rhs hPiii
Jn .s hxy i, speciesy (P), JIn (r).thxz i, speciesz (Q)
Jn .s hxy i, JIn (r).thxz i

! thxz i, speciesz (Q), s hwy i, speciesy (P), rhzw i

! thxz i, rhzw i, s hwy i

species(n · K) = K 7! n

w fresh

91

M. Miculan and I. Sambarino

The function species is used to initialise the abstract machine at the beginning of a simulation. If we aim
to simulate the execution of a process P ∈ Proc, the corresponding initial state (rooted in root) is
LPMroot , f n(J) ` J ⊕ (0, 0,
/ 0)
/

where J = species0,root
(P).
/

The reactions function is used to adjust the set of possible reactions dynamically.

5

Implementing the Stochastic Brane Calculus in COWGSAM

In this section, we show how the COW Generic Stochastic Abstract Machine can be used to implement
the Stochastic Brane Calculus, following the protocol described in the previous section.

5.1

Encoding of the Stochastic Brane Calculus

Syntax of species We define the species for the brane calculus, which in turn lead us to introduce
complexes and actions. Notice that (despite the deceiving syntax) species are not systems and complexes
are not membranes; nevertheless, actions are a subset of membranes.
I ::= Chxy i

(Species)

C ::= A1 , . . . , An
A ::= Jn .σ |

(Complexes)

JIn (τ).σ

| Kn .σ |

KIn .σ

| Gn (τ).σ

(Actions)

Node names can appear in the species in S and in reactions R.
The species function We can now provide the definition of the translation of a process P ∈ Bsys into a
species set. Basically, each compartment is assigned a different, fresh node name; therefore, the function
speciesE,x (P) is parametric in the set E of allocated node names, and the name x ∈ E to be used as the
location of the system P. The name x changes as we descend the compartment hierarchy.
In order to capture correctly the multiplicity of each species, we assume that systems are in normal
form. Basically, this form is a shorthand for products where n copies of the same system, i.e., P ◦ · · · ◦ P,
are represented as n · P.
Normal Systems :: Bnsys

Q ::= n1 · σ1 hQ1 i m . . . m nk · σk hQk i

where k ≥ 0 and for i 6= j : σi hQi i 6≡ σ j hQ j i

A system in normal form can be translated into a system just by unfolding the products. For Q ∈ Bnsys a
system in normal form, let dQe ∈ Bsys defined as follows:
n times

dke = k

z
}|
{
dn · σ hQ i ◦ Q e = σ hdQ0 ei ◦ · · · ◦ σ hdQ0 ei ◦dQ00 e
0

00

Proposition 5.1. For all P ∈ Bsys , there exists a system in normal form Q ∈ Bnsys such that dQe ≡ P.
As a consequence, we can give the definition of species on systems in normal form, as follows:
speciesE,x (k) , 0/
speciesE,x (n · σ hQ1 i ◦ Q2 ) , {s(σ )hxy i 7→ n} ∪ speciesE]{y},y (Q1 ) ∪ speciesE,x (Q2 )

with y ∈
/ E and f n(speciesE]{y},y (Q1 )) ∩ f n(speciesE,x (Q2 )) ⊆ {x}

92

Implementing Stochastics Brane Calculus in COWGSAM

The condition in the second case ensures that two different compartments are never given the same
name—any name clash has to be resolved by an α-conversion. The function s( ) converts a membrane
into a set of complexes:
s(0) , 0/

s(σ |σ 0 ) , s(σ ) ∪ s(σ 0 )

s(π.σ ) , {π.σ }

The reactions function The next step is to define the function reactions(I 7→ i, S), for I a species with
multiplicity i and S a population.
reactionsE (I1 7→ i, S) , unaryE (I1 ) ∪ binaryE (I1 , S)

unaryE (I1 ) , {({I1 7→ 1}, rn , id, {(s(σ ) ∪U10 )hxy i 7→ 1, s(ρ)hyw i 7→ 1}) |
I1 = U1 hxy i, U1 = {Gn (ρ).σ } ∪U10 , w ∈
/ E}


binaryE (I1 , S) , ({I1 7→ 1, I2 7→ 1}, rn , f , {(s(τ) ∪U10 ∪ s(σ ) ∪U20 )hxy i 7→ 1}) | I2 ∈ dom(S),
(I1 = U1 hxy i, I2 = U2 hyw i) ∨ (I2 = U1 hxy i, I1 = U2 hyw i),

U1 = {KIn .τ} ∪U10 ,U2 = {Kn .σ } ∪U20 , f = λ T.T [w := x] ]

({I1 7→ 1, I2 7→ 1}, rn , id,

{(s(τ) ∪U10 )hxy i 7→ 1, s(ρ)hyw i 7→ 1, (s(σ ) ∪U20 )hwz i 7→ 1}) | I2 ∈ dom(S),
(I1 = U1 hxy i, I2 = U2 hxz i) ∨ (I2 = U1 hxy i, I1 = U2 hxz i),
U1 = {JIn (ρ).τ} ∪U10 ,U2 = {Jn .σ } ∪U20 , w ∈
/E

In the case of Brane Calculus, the unary reactions are only those arising from pinocytosis, while binary
reactions arise from exocytosis and phagocytosis. In both cases, the multiplicity of each reactant is 1, so
the multiplicity of I1 is not relevant. Exocytosis merges two compartments; this is reflected by the fact
that the “rearranging” function f substitutes every occurrence of the name w in T with x. On the other
hand, pinocytosis and phagocytosis create new compartments; to represent the new structure, we choose
a fresh name w representing the new intermediate nesting level, and reconnect the various compartments
accordingly. Therefore, for any reaction (S1 , r, f , S2 ) ∈ reactionsE (I1 7→ i, S), f n(S2 ) \ E is either 0/ (in
the case of exocytosis) or the singleton {w}.

5.2

Adequacy results

Before proving the correctness of our implementation, we have to define how to translate a species set
back to a system of the brane calculus.
Let S be a non empty species set. A root name of S, denoted by root(S), is a name x such that
S(Chxy i) > 0 for some C, y, and for all z,C0 : C0 hzx i ∈
/ dom(S). The next result states that root( ) is well
defined on the species sets we encounter during a simulation.
Lemma 5.2. For all P ∈ Bsys :

1. if P 6= k: root(LPMx ) = x.
a,O

2. if LPMx −−→ E ` (t, S, R) and S 6= 0,
/ then root(S) = x.

Proof. (1.) is trivial by definition. (2.) It is enough to check that the reaction rules introduced by
reactions( ) do not change the name of the root, nor introduce new ones.

93

M. Miculan and I. Sambarino

We can now define a function J K which maps complexes to membranes, and a function J Kx mapping
species sets to systems; the latter is parametric in the name x of the root of the system.
JA1 , . . . , An K , A1 | . . . |An

JSKx ,

∏

Chxy i∈dom(S)

S(Chxy i) · (JCKhJSKy i)

where the notation n · P is a shorthand for P m . . . m P, n times.
Lemma 5.3. For all P ∈ Bsys :
1. JLPMx Kx ≡ P.
a,O

2. if LPMx −−→ E ` (t, S, R) then JSKx is well defined.

Proof. (1.) is easy. (2.) It is enough to check that the reaction rules introduced by reactions( ) do not
introduce loops (i.e., the order among names is well founded).
We can now state and prove the main results of this section.
a,O

Proposition 5.4 (Soundness). For all P ∈ Bsys , if LPMx −−→ E ` (t, S, R) then there exists µ such that
P→
− µ and µid ([JSKx ]) = a.

Proof. The proof is by cases on which reaction rule O is selected by the function next. By additivity of
measures, we can restrict ourselves to when the whole process P is the redex of the reduction. Let us see
here the case when P is the redex of a (red-pino) (another is in Appendix B).
Let P = Gn (ρ).σ |σ0 hP0 i and let us assume that σ0 does not exhibit a Gn action. Then, the translation
of P is LPMx = E ` species0,x
/ 0),
/ where
/ (P) ⊕ (0, 0,
x
0
species0,x
/ (P) = {s(Gn (ρ).σ |σ0 )hy i 7→ 1} ∪ species{y},x (P )

= {(s(Gn (ρ).σ ) ∪ s(σ0 ))hxy i 7→ 1} ∪ species{y},x (P0 )
= {({Gn (ρ).σ } ∪ s(σ0 ))hxy i 7→ 1} ∪ species{y},x (P0 )

Let I1 = ({Gn (ρ).σ } ∪ s(σ0 ))hxy i and JP0 = species{y},x (P0 ); then LPMx = I1 7→ 1 ⊕ JP0 ⊕ (0, 0,
/ 0)
/ = JP0 ⊕
0
0
(0, S , R ) where
S0 = {I1 7→ 1}

R0 = init(L, (0, S0 , 0))
/ = {OL 7→ (t1 , a1 )}

OL = ({I1 7→ 1}, rn , id, {(s(σ ) ∪ s(σ0 ))hxy i 7→ 1, s(ρ)hyw i 7→ 1})

L = reactions(I1 7→ 1, 0)
/ = unary(I1 ) = {({I1 7→ 1}, rn , id, {(s(σ ) ∪ s(σ0 ))hxy i 7→ 1, s(ρ)hyw i 7→ 1})}
a,O

Now, the reaction O in LPMx −−→ T is OL (otherwise it would involve P0 , not the pino of the whole P).
a,O

This means that LPMx −−→ T is derived by means of an application of the (Reaction rule) as follows, where
S1 = {I1 7→ 1} and S2 = {(s(σ ) ∪ s(σ0 ))hxy i 7→ 1, s(ρ)hyw i 7→ 1}:
((S1 , rn , id, S2 ), a1 ,t1 ) = next(0, S0 , R0 ) (E 0 ` (0, S00 , R00 )) = cow(E ` (0, S0 , R0 ), S1 )
a1 ,(S1 ,rn ,id,S2 )

E ` {JP0 } ⊕ (0, S0 , R0 ) −−−−−−−−→ norm(E ∪ f n(S2 ) ` ({JP0 } ⊕ S2 ⊕ ((t1 , S00 , R00 )
where S00 = S0 and R00 = R0 .

S1 )))

94

Implementing Stochastics Brane Calculus in COWGSAM
{JP0 } ⊕ S2 ⊕ ((t1 , S00 , R00 ) S1 ) =
= JP0 ⊕ {(s(σ ) ∪ s(σ0 ))hxy i 7→ 1, s(ρ)hyw i 7→ 1} ⊕ ((t1 , S0 , R0 ) {I1 7→ 1})
= JP0 ⊕ {(s(σ ) ∪ s(σ0 ))hxy i 7→ 1, s(ρ)hyw i 7→ 1} ⊕ (t1 , {I1 7→ 0}, {OL 7→ (t2 , a2 )})
= JP0 ⊕ J 0 ⊕ (t1 , {I1 7→ 0}, {OL 7→ (t2 , a2 )})
Now let us define Q as Q = σ |σ0 hρhki m P0 i, then LQMx = species0,x
/ 0)
/ where
/ (Q) ⊕ (0, 0,
0
species0,x
/ (Q) = species0,x
/ (σ |σ0 hρhki m P i)

= {s(σ |σ0 )hxy i 7→ 1} ∪ species{x},y (ρhki m P0 )

= {(s(σ ) ∪ s(σ0 ))hxy i 7→ 1} ∪ species{x},y (ρhki) ∪ species{x},y (P0 )
= {(s(σ ) ∪ s(σ0 ))hxy i 7→ 1} ∪ {s(ρ)hyw i 7→ 1} ∪ 0/ ∪ JP0 = J 0 ∪ JP0

and hence Q ≡ JJP0 ⊕ J 0 ⊕ (t1 , {I1 7→ 0}, {OL 7→ (t2 , a2 )})Kx . It remains to prove that rn = µid ([Q]). Let
us notice that the derivation of P → µ is actually as follows:
ρ

(par)
(loc)
G (ρ).σ |σ0

where µ = µ 0 @P0n

σ0 →
− µ 00
ρ
Gn (ρ).σ |σ0 →
− [Gn ]σ Gn (ρ).σ σ0 µ 00
Gn (ρ).σ |σ0 hP0 i →
− µ
Gn (ρ).σ →
− [Gn ]σ

ρ

P0 →
− µ0

([Gn ]σ Gn (ρ).σ σ0 µ 00 ). Then:
G (ρ).σ |σ0

µid ([σ |σ0 hρhki m P0 i]) = (µ 0 @P0n
=

=
=
=

ρ

00
0
Gn (ρ).σ σ0 µ ))id ([σ |σ0 hρhki m P i])
ρ
µid0 ([σ |σ0 hρhki m P0 i]) + ([Gn ]σ Gn (ρ).σ σ0 µ 00 )Gn ([σ |σ0 ] × [ρ])
ρ
([Gn ]σ Gn (ρ).σ σ0 µ 00 )Gn ([σ |σ0 ] × [ρ])
ρ
([Gn ]σ )Gn ([σ |σ0 ] × [ρ]) + µG00 n ([σ |σ0 ] × [ρ])
rn + µG00 n ([σ |σ0 ] × [ρ]) = rn

([Gn ]σ

where the last equivalence holds because µG00 n ([σ |σ0 ] × [ρ]) = 0 because we assumed that the reaction
does not involve σ0 .
Proposition 5.5 (Progress). For all processes P, Q, if P } Q then there exists a reaction O and a term
a,O

T such that LPMx −−→ T and Q ≡ JT Kx .

Proof. By induction on the derivation of P}Q. Let us see the case of (red-pin), the others being similar.
Let P = G(ρ).σ |σ0 hP0 i and Q = σ |σ0 hρhki m P0 i. Then,
LPMx = species{y},x (P0 ) ⊕ (0, {({Gn (ρ).σ } ∪ s(σ0 ))hxy i 7→ 1}, {OL 7→ (t1 , a1 )})

where OL = ({({Gn (ρ).σ } ∪ s(σ0 ))hxy i 7→ 1}, rn , id, {(s(σ ) ∪ s(σ0 ))hxy i 7→ 1, s(ρ)hyw i 7→ 1}). Then, by
the (Reaction rule) we can take T = species{x},y (P0 ) ⊕ (t1 , {(s(σ ) ∪ s(σ0 ))hxy i 7→ 1, s(ρ)hyw i 7→ 1}, 0).
/ It
is easy to check that Q ≡ JT Kx .

Proposition 5.6 (Completeness). For all processes P, Q, if P →
− µ and µid ([Q]) > 0 then for all node
a,O

name x, there exists a reaction O and a term T such that LPMx −−→ T , Q ≡ JT Kx and a = µid ([Q]).

Proof. If P →
− µ and µid ([Q]) > 0 then P } Q by Prop. 3.4. By Prop. 5.5, we have that for some a, O, T ,
a,O

LPMx −−→ T and Q ≡ JT Kx . But then a = µid ([Q]) by soundness (Prop. 5.4).

M. Miculan and I. Sambarino

5.3

95

Example

We conclude this section with an example. Let P = 10000·Jn .Km hJk hiim100·((JIn (KIm ) | KI )hJk hii).
Then, its reductions in the Brane Calculus are as follows:
10000 · Jn .Km hJk hii m 100 · ((JIn (KIm ) | KI )hJk hii)

} 9999 · Jn .Km hJk hii m 99 · ((JIn (KIm ) | KI )hJk hii) m KI hKIm hKm hJk hiii m Jk hii
} 9999 · Jn .Km hJk hii m 99 · ((JIn (KIm ) | KI )hJIk ()hii) m KI h2 · Jk hi m khkii

The translation of P is LPMx = E ` species0,x
/ 0),
/ where
/ (P) ⊕ (0, 0,

I
I
I
species0,x
/ (P) = species{x} (10000 · Jn .Km hJk hii ◦ 100 · (Jn (Km ) | K )hJk hii)

= {s(Jn .Km )hxy i 7→ 10000} ∪ species{x,y},y (Jk hi) ∪ {s(JIn (KIm ) | KI )hxz i 7→ 100}
∪ species{x,z},z (Jk hi)

= {{Jn .Km }hxy i 7→ 10000, {JIn (KIm ), KI }hxz i 7→ 100, {Jk }hyw i 7→ 1,
{Jk }hzv i 7→ 1}

Let I1 = {Jn .Km }hxy i, I2 = {JIn (KIm ), KI }hxz i, I3 = {Jk }hyw i, I4 = {Jk }hzv i, rn = 10s−1 , rk = 5s−1 and
rm = 5s−1 ; then
LPMx = I1 7→ 10000 ⊕ I2 7→ 100 ⊕ I3 7→ 1 ⊕ I4 7→ 1 ⊕ (0, 0,
/ 0)
/
= I2 7→ 100 ⊕ I3 7→ 1 ⊕ I4 7→ 1 ⊕ (0, S1 , R1 )
= I3 7→ 1 ⊕ I4 7→ 1 ⊕ (0, S2 , R2 )

= I4 7→ 1 ⊕ (0, S3 , R3 ) = (0, S4 , R4 )

where
L1 = reactions(I1 7→ 10000, 0)
/ = 0/
S1 = S{I1 7→ 10000}

R1 = R ∪ init(L1 , (0, S1 , 0))
/ = 0/
L2 = reactions(I2 7→ 100, S1 )

= ({I2 7→ 1, I1 7→ 1}, 10, id, {{KI }hxz i 7→ 1, {KIm )}hzw i 7→ 1, {Km }hwy i 7→ 1})

S2 = S1 {I2 7→ 100} = {I1 7→ 10000, I2 7→ 100}

R2 = R1 ∪ init(L2 , (0, S2 , R1 )) = {OL2 7→ (t1 , a1 )}
L3 = reactions(I3 7→ 1, S2 ) = 0/

S3 = S2 {I3 7→ 1} = {I1 7→ 10000, I2 7→ 100, I3 7→ 1}

R3 = R2 ∪ init(L3 , (0, S3 , R2 )) = R2
L4 = reactions(I4 7→ 1, S3 ) = 0/

S4 = S3 {I4 7→ 1} = {I1 7→ 10000, I2 7→ 100, I3 7→ 1, I4 7→ 1}

R4 = R3 ∪ init(L4 , (0, S4 , R3 )) = R3

96

Implementing Stochastics Brane Calculus in COWGSAM

with OL2 = ({I2 7→ 1, I1 7→ 1}, 10, id, {{KI }hxz i 7→ 1, {KIm }hzw i 7→ 1, {Km }hwy i 7→ 1})
a1 = propensity(OL2 , S2 ) = 10000000
t1 = 0 + delay(10, 10000000).
a,O

a,O

Now, the reaction O in LPMx −−→ T is OL . This means that LPMx −−→ T is derived by means of an
application of the (Reaction rule) as follows, where E = x, y, z, w, S1 = {I2 7→ 1, I1 7→ 1}, and S2 =
{{KI }hxz i 7→ 1, {KIm }hzw i 7→ 1, {Km }hwy i 7→ 1}:
((S1 , 10, id, S2 ), 10,t1 ) = next(0, S4 , R4 ) (E 0 ` (0, S5 , R5 )) = cow(E ` (0, S4 , R4 ), S1 )
10,(S1 ,10,id,S2 )

E ` (0, S4 , R4 ) −−−−−−−−→ norm(E 0 ∪ f n(S2 ) ` (S2 ⊕ ((t1 , S5 , R5 )

S1 )))

where
S5 = {I1 7→ 1, {Jn .Km }hxy0 i 7→ 9999, I2 7→ 1, {JIn (KIm ), KI }hxz0 i 7→ 99, I3 7→ 1,
0

I4 7→ 1, {Jk }hyw0 i 7→ 1, {Jk }hzv0 i 7→ 1}
0

R5 = {OL2 7→ (t1 , a1 ), O1 7→ (t2 , a2 ), O2 7→ (t3 , a3 ), O3 7→ (t4 , a4 )}
with
00

O1 = ({{Jn .Km }hxy0 i 7→ 1, I2 7→ 1}, 10, id, {{KI }hxz i 7→ 1, {KIm }hzw00 i 7→ 1, {Km }hwy0 i 7→ 1})
0

000

O2 = ({JIn (KIm ), KI }hxz0 i 7→ 1, I1 7→ 1}, 10, id, {{KI }hxz0 i 7→ 1, {KIm }hzw000 i 7→ 1, {Km }hwy i 7→ 1})
0

00

O3 = ({{Jn .Km }hxy0 i 7→ 1, {JIn (KIm ), KI }hxz0 i 7→ 1}, 10, id, {{JI }hxz0 i 7→ 1, {KIm }hzz00 i 7→ 1, {Km }hzy0 i})
a2 = propensity(O1 , S5 ) = 99990
a3 = propensity(O2 , S5 ) = 990
a4 = propensity(O3 , S5 ) = 989901
t2 = 0 + delay(10, 99990)
t3 = 0 + delay(10, 990)
t4 = 0 + delay(10, 989901)
We can now compute the multiset of the new machine state:
S2 ⊕ ((t1 , S5 , R5 )

S1 ) = S2 ⊕ ((t1 , S5 , R5 )

= S2 ⊕ ((t1 , S6 , R6 })
= {{K

I

}hxz i

= (t1 , S8 , R8 )

7→

{I2 7→ 1, I1 7→ 1})

{I1 7→ 1})
I
1, {Km }hzw i 7→ 1, {Km }hwy i

7→ 1} ⊕ (t1 , S7 , R7 )

with
a02 = propensity(O1 , S6 ) = 0

S8 = S6 ∪ {{KI }hxz i 7→ 1, {KIm }hzw i 7→ 1, {Km }hwy i 7→ 1}

t20 = t1 + (a2 /a02 )(t2 − t1 )

R6 = {O1 7→ (t20 , a02 )}

t30 = t1 + (a3 /a5 )(t3 − t1 )

R8 = R7 ∪ {{{KIm }hzw i 7→ 1, {Km }hwy i 7→ 1}, 5, f , {0h
/ zw i} 7→ 1}

a03 = propensity(O2 , S7 ) = 0
S6 = S5 \{I2 7→ 1}

S7 = S6 \{I1 7→ 1}

R7 = R6 ∪ {O2 7→ (t30 , a03 )}

f = λ T.T [y := z]

M. Miculan and I. Sambarino

6

97

Conclusions

In this paper, we have presented an abstract machine for the Stochastic Brane Calculus. Instead of defining an ad hoc machine, we have adopted the generic abstract machine for stochastic calculi (GSAM) recently introduced by Lakin, Paulevé and Phillips. According to the encoding technique we have adopted,
membranes are flattened into a set of species, where the hierarchical structure is represented by means
of names. In order to keep track of these names, and for dealing efficiently with multiple copies of the
same species, we have introduced a new generic abstract machine, called COWGSAM, which extends
the GSAM with a name environment and a copy-on-write optimization strategy. We have proved that
the implementation of the Stochastic Brane Calculus in COWGSAM is adequate with respect to the
stochastic structural operational semantics of the calculus given in [2].
We think that COWGSAM can be used for implementing other stochastic calculi dealing with nested
structures, also beyond the models for systems biology. In particular, it is interesting to apply this approach to Stochastic Bigraphs [9], a general meta-model well-suited for representing a range of stochastic systems with compartments; in this way we would obtain a General Stochastic Bigraphical Machine,
which could be instantiated to any given stochastic bigraphic reactive system. However, such a machine
would not scale well, as in general the COW strategy may be not very useful; thus, we can restrict our
attentions to smaller subsets of BRSs, specifically designed to some application domain. For biological
applications, the bigraphic reactive systems considered in [1, 7] might be a more reasonable target.
Another interesting question is about the expressive power of GSAM and COWGSAM. We think
that GSAM correspond to stochastic (multiset) Petri nets, but COWGSAM could go further thanks to
the possibility of creating unlimited new names during execution. Further work include comparison with
other stochastic simulation tools dealing with compartments, like BioPEPA [6].
Acknowledgment Work funded by MIUR PRIN project “SisteR”, prot. 20088HXMYN.

References
[1] Giorgio Bacci, Davide Grohmann & Marino Miculan (2009): Bigraphical models for protein and membrane
interactions. In: Gabriel Ciobanu, editor: Proc. MeCBIC’09, EPTCS 11, pp. 3–18. doi:10.4204/EPTCS.11.1.
[2] Giorgio Bacci & Marino Miculan (2012): Measurable Stochastics for Brane Calculus. Theoretical Computer
Science 431, pp. 117–136. doi:10.1016/j.tcs.2011.12.055.
[3] Roberto Barbuti, Andrea Maggiolo-Schettini, Paolo Milazzo & Angelo Troina (2007): The Calculus of Looping Sequences for Modeling Biological Membranes. In: George Eleftherakis, Petros Kefalas, Gheorghe Paun,
Grzegorz Rozenberg & Arto Salomaa, editors: Workshop on Membrane Computing, Lecture Notes in Computer Science 4860, Springer, pp. 54–76. doi:10.1007/978-3-540-77312-2 4.
[4] Luca Cardelli (2004): Brane Calculi. In: Vincent Danos & Vincent Schächter, editors: Proc. CMSB, Lecture
Notes in Computer Science 3082, Springer, pp. 257–278. doi:10.1007/978-3-540-25974-9 24.
[5] Luca Cardelli & Radu Mardare (2010): The Measurable Space of Stochastic Processes. In: Proc. QEST,
IEEE Computer Society, pp. 171–180. doi:10.1109/QEST.2010.30.
[6] Federica Ciocchetta & Maria Luisa Guerriero (2009): Modelling Biological Compartments in Bio-PEPA.
Electronic Notes in Theoretical Computer Science 227, pp. 77–95. doi:10.1016/j.entcs.2008.12.105.
[7] Troels Christoffer Damgaard, Espen Højsgaard & Jean Krivine (2012): Formal Cellular Machinery. Electronic Notes in Theoretical Computer Science 284, pp. 55–74. doi:10.1016/j.entcs.2012.05.015.
[8] Holger Hermanns (2002): Interactive Markov Chains: The Quest for Quantified Quality, Lecture Notes in
Computer Science 2428. Springer. doi:10.1007/3-540-45804-2.

98

Implementing Stochastics Brane Calculus in COWGSAM

[9] Jean Krivine, Robin Milner & Angelo Troina (2008): Stochastic Bigraphs. In: Proc. 24th MFPS, ENTCS
218, pp. 73–96. doi:10.1016/j.entcs.2008.10.006.
[10] Matthew R. Lakin, Loı̈c Paulevé & Andrew Phillips (2012): Stochastic Simulation of Multiple Process Calculi for Biology. Theoretical Computer Science 431, pp. 181–206. doi:10.1016/j.tcs.2011.12.057.
[11] Cosimo Laneve & Fabien Tarissan (2008): A simple calculus for proteins and cells. Theoretical Computer
Science 404(1-2), pp. 127–141. doi:10.1016/j.tcs.2008.04.011.
[12] Prakash Panangaden (2009): Labelled Markov Processes. Imperial College Press, London, U.K.
[13] Loı̈c Paulevé, Simon Youssef, Matthew R. Lakin & Andrew Phillips (2010): A generic abstract machine for stochastic process calculi. In: Paola Quaglia, editor: Proc. CMSB, ACM, pp. 43–54.
doi:10.1145/1839764.1839771.
[14] Aviv Regev, Ekaterina M. Panina, William Silverman, Luca Cardelli & Ehud Y. Shapiro (2004): BioAmbients: an abstraction for biological compartments. Theoretical Computer Science 325(1), pp. 141–167.
doi:10.1016/j.tcs.2004.03.061.

A

Some measure theory

Given a set M, a family Σ of subsets of M is called a σ -algebra if it contains M and is closed under
complements and (infinite) countable unions:
1. M ∈ Σ;

2. A ∈ Σ implies Ac ∈ Σ, where Ac = M \ A;

3. {Ai }i∈N ⊂ Σ implies

i∈N Ai

S

∈ Σ.

Since M ∈ Σ and M c = 0,
/ 0/ ∈ Σ, hence Σ is nonempty by definition. A σ -algebra is closed under countable
set-theoretic operations: is closed under finite unions (A, B ∈ Σ implies A ∪ B = A ∪ B ∪ 0/ ∪ 0/ ∪ · · · ∈ Σ),
countable intersections (by DeMorgan’s law A ∩ B = (Ac ∪ Bc )c in its finite and inifite version), and
countable subtractions (A, B ∈ Σ implies A \ B = A ∩ Bc ∈ Σ).
Definition A.1 (Measurable Space). Given a set M and a σ -algebra on M, the tuple (M, Σ) is called a
measurable space, the elements of Σ measurable sets, and M the support-set.

A set Ω ⊆ 2M is a generator for the σ -algebra Σ on M if Σ is the closure of Ω under complement
and countable union; we write σ (Ω) = Σ and say that Σ is generated by Ω. Note that the σ -algebra generated by a Ω is also the smallest σ -algebra containing Ω, that is, the intersection of all σ -algebras that
contain Ω. In particular it holds that a completely arbitrary intersection of σ -algebras is a σ -algebra. A
σ -algebra generated by Ω, denoted by σ (Ω), is minimal in the sense that if Ω ⊂ Σ and Σ is a σ -algebra,
then σ (Ω) ⊂ Σ. If Ω is a σ -algebra then obviously σ (Ω) = Ω; if Ω is empty or Ω = {0},
/ or Ω = {M},
then σ (Ω) = {0,
/ M}; if Ω ⊂ Σ and Σ is a σ -algebra, then σ (Ω) ⊂ Σ. A generator Ω for Σ is a base for Σ
if it has disjoin elements. Note that if Ω is a base for Σ, all measurable sets in Σ can be decomposed into
countable unions of elements in Ω.
+
A measure on a measurable space (M, Σ) is a function µ : Σ → R+
∞ , where R∞ denotes the extended
positive real line, such that
1. µ(0)
/ = 0;
2. for any disjoint sequence {Ni }i∈I ⊆ Σ with I ⊆ N, it holds
µ(

i∈I Ni )

S

= ∑i∈I µ(Ni ) .

99

M. Miculan and I. Sambarino

The triple (M, Σ, µ) is called a measure space. A measure space (M, Σ, µ) is called finite if µ(M) is a
finite real number; it is called σ -finite if M can be decomposed into a countable union of measurable sets
S
of finite measure, that is, M = i∈I Ni , for some I ⊆ N and µ(Ni ) ∈ R+ for each i ∈ I. A set in a measure
space has σ -finite measure if it is a countable union of sets with finite measure. Specifying a measure
includes specifying its domain. If µ is a measure on a measurable space (M, Σ) and Σ0 is a σ -algebra
contained in Σ, then the restriction µ 0 of µ to Σ0 is also a measure, and in particular a measure on (M 0 , Σ0 ),
for some M 0 ⊆ M such that Σ0 is a σ -algebra on M 0 .
Given two measurable spaces and measures on them, one can obtain the product measurable space
and the product measure on that space. Let (M1 , Σ1 ) and (M2 , Σ2 ) be measurable spaces, and µ1 and
µ2 be measures on these spaces. Denote by Σ1 ⊗ Σ2 the σ -algebra on the cartesian product M1 × M2
generated by subsets of the form B1 × B2 , said rectangles, where B1 ∈ Σ1 and B2 ∈ Σ2 . The product
measure µ1 ⊗ µ2 is defined to be the unique measure on the measurable space (M1 × M2 , Σ1 ⊗ Σ2 ) such
that, for all B1 ∈ Σ1 and B2 ∈ Σ2
(µ1 ⊗ µ2 )(B1 × B2 ) = µ1 (B1 ) · µ2 (B2 )
The existence of this measure is guaranteed by the Hahn-Kolmogorov theorem. The uniqueness of the
product measure is guaranteed only in the case that both (M1 , Σ1 , µ1 ) and (M2 , Σ2 , µ2 ) are σ -finite.
Let ∆(M, Σ) be the family of measures on (M, Σ). It can be organized as a measurable space by
considering the σ -algebra generated by the sets {µ ∈ ∆(M, Σ) : µ(S) ≥ r}, for arbitrary S ∈ Σ and r > 0.
Given two measurable spaces (M, Σ) and (N, Θ) a mapping f : M → N is measurable if for any
T ∈ Θ, f −1 (T ) ∈ Σ. Measurable functions are closed under composition: given f : M → N and g : N → O
measurable functions then g ◦ f : M → O is also measurable.

B

Proof of Prop. 5.4

Let P = KIn .τ|τ0 hKn .σ |σ0 hP0 i m P00 i; then, LPMx = E ` species0,x
/ 0),
/ where
/ (P) ⊕ (0, 0,
I
x
0
00
species0,x
/ (P) = {s(Kn .τ|τ0 )hy i} ∪ species{x},y (Kn .σ |σ0 hP i m P )

= {({KIn .τ} ∪ s(τ0 ))hxy i 7→ 1} ∪ species{x},y (Kn .σ |σ0 hP0 i) ∪ species{x},y (P00 )

= {({KIn .τ} ∪ s(τ0 ))hxy i 7→ 1} ∪ {s(Kn .σ |σ0 )hyw i 7→ 1} ∪ species{x,y},w (P0 ) ∪ species{x},y (P00 )

= {({KIn .τ} ∪ s(τ0 ))hxy i 7→ 1} ∪ {(s(Kn .σ ) ∪ s(σ0 ))hyw i 7→ 1} ∪ species{x,y},w (P0 ) ∪ species{x},y (P00 )

= {({KIn .τ} ∪ s(τ0 ))hxy i 7→ 1} ∪ {({Kn .σ } ∪ s(σ0 ))hyw i 7→ 1} ∪ species{x,y},w (P0 ) ∪ species{x},y (P00 )

Let I1 = ({KIn .τ} ∪ s(τ0 ))hxy i, I2 = ({Kn .σ } ∪ s(σ0 ))hyw i, JP0 = species{x,y},w (P0 ), JP00 = species{x},y (P00 );
then LPMx = I1 7→ 1 ⊕I2 7→ 1 ⊕JP0 ⊕JP00 ⊕(0, 0,
/ 0)
/ = I2 7→ 1 ⊕JP0 ⊕JP00 ⊕(0, S0 , R0 ) = JP0 ⊕JP00 ⊕(0, S00 , R00 )
where
L0 = reactions(I1 7→ 1, 0)
/ = 0/
S0 = S{I1 7→ 1}

R0 = init(L0 , (0, S0 , 0))
/ = 0/
L00 = reactions(I2 7→ 1, S0 )

= ({I2 7→ 1, I1 7→ 1}, rn , f , {(s(τ) ∪ s(τ0 ) ∪ s(σ ) ∪ s(σ0 ))hxy i 7→ 1})

S00 = S0 {I2 7→ 1} = {I1 7→ 1, I2 7→ 1}

R00 = init(L00 , (0, S00 , R0 )) = {OL 7→ (t1 , a1 )}

100

Implementing Stochastics Brane Calculus in COWGSAM

with OL = ({I1 7→ 1, I2 7→ 1}, rn , f , {(s(τ) ∪ s(τ0 ) ∪ s(σ ) ∪ s(σ0 ))hxy i 7→ 1}). Now, the reaction O in
F,O

F,O

LPMx −−→ T is OL . This means that LPMx −−→ T is derived by means of an application of the (Reaction rule)
as follows, where S1 = {I1 7→ 1, I2 7→ 1}, S2 = {(s(τ) ∪ s(τ0 ) ∪ s(σ ) ∪ s(σ0 ))hxy i 7→ 1} and f = T [w := x]:
((S1 , rn , f , S2 ), a1 ,t1 ) = next(0, S00 , R00 ) (E 0 ` (0, S000 , R000 )) = cow(E ` (0, S00 , R00 ), S1 )
a1 ,(S1 ,rn , f ,S2 )

where

S000

=

E ` (0, S00 , R00 ) −−−−−−−−→ norm(E 0 ∪ f n(S2 ) ` ( f (S2 ⊕ ((t1 , S000 , R000 )

S00

and

R000

( f (S2 ⊕ ((t1 , S000 , R000 )

=

S1 ))))

R00 .

S1 )))) =

= f ((s(τ) ∪ s(τ0 ) ∪ s(σ ) ∪ s(σ0 )hxy i) ⊕ ((t1 , S000 , R000 )

{I1 7→ 1, I2 7→ 1})

= f ((s(τ) ∪ s(τ0 ) ∪ s(σ ) ∪ s(σ0 )hxy i) ⊕ ((t1 , {I1 7→ 0, I2 7→ 1}, {OL 7→ (t2 , a2 )})
=

=
=

f ((s(τ) ∪ s(τ0 ) ∪ s(σ ) ∪ s(σ0 )hxy i) ⊕ (t1 , {I1 7→ 0, I2 7→ 0}, {OL 7→ (t3 , a3 )})
(s(τ) ∪ s(τ0 ) ∪ s(σ ) ∪ s(σ0 )hxy i) ⊕ (t1 , {I1 7→ 0, I2 7→ 0}, {OL 7→ (t3 , a3 )})
J 0 ⊕ (t1 , {I1 7→ 0, I2 7→ 0}, {OL 7→ (t3 , a3 )})

{I2 7→ 1})

Now let us define Q as Q = σ |σ0 |τ|τ0 hP00 i m P0 , then LQMx = species0,x
/ 0)
/ where
/ (Q) ⊕ (0, 0,
00
0
species0,x
/ (Q) = species0,x
/ (σ |σ0 |τ|τ0 hP i m P )

00
0
= species0,x
/ (σ |σ0 |τ|τ0 hP i) ∪ species0,x
/ (P )

0
= {s(σ |σ0 |τ|τ0 )hxy i) 7→ 1} ∪ species{x},y (P00 ) ∪ species0,x
/ (P )

= {s(τ) ∪ s(τ0 ) ∪ s(σ ) ∪ s(σ0 )hxy i 7→ 1} ∪ JP00 ∪ JP0 = J 0 ∪ JP00 ∪ JP0

And hence Q ≡ JJP0 ⊕ JP00 ⊕ J 0 ⊕ (t1 , {I1 7→ 0, I2 7→ 0}, {OL 7→ (t3 , a3 )})Kx . It remains to prove that rn =
µid ([Q]). Let us notice that the derivation of P →
− µ is actually as follows:
Kn .σ →
− [Kn ]σ
σ0 →
− µ 00
(par)
Kn .σ |σ0 →
− [Kn ]σ Kn .σ σ0 µ 00

P0 →
− µ0

(loc)
0 − µ 0 @Kn .σ |σ0 ([K ]
00
P00 →
− µ 000
n σ Kn .σ σ0 µ )
P0
KIn .τ →
− [KIn ]τ
τ0 →
− µ (comp) K.σ |σ0 hP i →
(par)
K .σ |σ
KIn .τ|τ0 →
− [KIn ]τ KIn .τ τ0 µ
Kn .σ |σ0 hP0 i m P00 →
− (µ 0 @P0n 0 ([Kn ]σ Kn .σ σ0 µ 00 ))Kn .σ |σ0 hP0 i ⊗P00 µ 000
(loc)
KIn .τ|τ0 hKn .σ |σ0 hP0 i m P00 i →
− ν
K .σ |σ0

where ν = ((µ 0 @P0n
K .σ |σ0

µ1 = (µ 0 @P0n

KI .τ|τ

([Kn ]σ Kn .σ σ0 µ 00 ))Kn .σ |σ0 hP0 i ⊗S µ 000 )@Knn .σ |σ00 hP0 imP00 ([KIn ]τ KIn .τ τ0 µ),

([Kn ]σ Kn .σ σ0 µ 00 ))Kn .σ |σ0 hP0 i ⊗P00 µ 000 and µ2 = [KIn ]τ KIn .τ τ0 µ. Then:
KI .τ|τ

νid ([σ |σ0 |τ|τ0 hP00 i m P0 ]) = (µ1 @Knn .σ |σ00 hP0 imP00 µ2 )id ([σ |σ0 |τ|τ0 hP00 i m P0 ])
= µ1 id ([σ |σ0 |τ|τ0 hP00 i m P0 ]) + µ1 exn ([σ |σ0 ] × [P0 ] × [P00 ]) · µ2 KIn ([τ|τ0 ])/rn

= µ1 exn ([σ |σ0 ] × [P0 ] × [P00 ]) · µ2 KIn ([τ|τ0 ])/rn
K .σ |σ0

= ((µ 0 @P0n
= (([Kn ]σ
=
=

000
([Kn ]σ Kn .σ σ0 µ 00 ))exn ([σ |σ0 ] × [P0 ] × [P00 ]) + µex
([σ |σ0 ] × [P0 ] × [P00 ])) · µ2 KIn ([τ|τ0 ])/rn
n

Kn .σ σ0

000
µ 00 )Kn ([σ |σ0 ]) + µex
([σ |σ0 ] × [P0 ] × [P00 ])) · ([KIn ]τ KIn .τ τ0 µ)KI ([τ|τ0 ])/rn
n
n

000
(([Kn ]σ )Kn ([σ |σ0 ]) + µK00 n ([σ |σ0 ]) + µex
([σ |σ0 ] × [P0 ] × [P00 ])) · (([KIn ]τ )KIn ([τ|τ0 ]) + µKIn ([τ|τ0 ]))/rn
n
000
(rn + µK00 n ([σ |σ0 ]) + µex
([σ |σ0 ] × [P0 ] × [P00 ])) · (rn + µKIn ([τ|τ0 ]))/rn = rn
n

000 ([σ |σ ] × [P0 ] × [P00 ]) = µ I ([τ|τ ]) = 0
where the last equivalence holds because µK00 n ([σ |σ0 ]) = µex
0
Kn
0
n
because we assumed that the reaction does not involve either σ0 nor τ0 .

