Exact rule-based stochastic simulations for the system with unlimited
number of molecular species
Anton V. Bernatskiy1, a) and Elizaveta A. Guseva2, b)
1)

Center for Complex Systems, University of Vermont, 210 Colchester Ave, Burlington, VT,
05405
2)
Laufer Center for Physical and Quantitative Biology, Stony Brook University, Stony Brook,
NY 17794-5252

arXiv:1609.06403v2 [] 4 Oct 2016

(Dated: 5 October 2016)

We introduce expandable partial propensity direct method (EPDM) – a new exact stochastic simulation
algorithm suitable for systems involving many interacting molecular species. The algorithm is especially
efficient for sparsely populated systems, where the number of species that may potentially be generated is
much greater than the number of species actually present in the system at any given time. The number of
operations per reaction scales linearly with the number of species, but only those which have one or more
molecules. To achieve this kind of performance we are employing a data structure which allows to add and
remove species and their interactions on the fly. When a new specie is added, its interactions with every
other specie are generated dynamically by a set of user-defined rules. By removing the records involving the
species with zero molecules, we keep the number of species as low as possible. This enables simulations of
systems for which listing all species is not practical. The algorithm is based on partial propensities direct
method (PDM) by Ramaswamy et al. for sampling trajectories of the chemical master equation.
Frequently used abbreviations and symbols
SSA – Stochastic Simulation Algorithm
DM – Direct Method17
PDM – Partial propensity Direct Method34
SPDM – Sorting Partial propensity Direct Method34
EPDM – Expandable Partial propensity Direct Method
N – number of molecular or agent species
M – number of all possible reactions or agent interactions
α – maximum number of possible reactions or agent interactions between any pair of molecular or agent species

I.

INTRODUCTION

Mathematical modeling of chemical reactions is an important part of systems biology research.
The traditional approach to modeling chemical systems is based on the law of mass of action equations
which assumes reactions to be macroscopic, continuous
and deterministic. It doesn’t account for the discreteness
of the reactions or temporal and spatial fluctuations, describing only the average properties instead. This makes
it ill-suited for modeling nonlinear systems or systems
with a small number of participating molecules, such as
living cells.
In contrast, stochastic simulation algorithms (SSAs) do
account for discreteness and inevitable randomness of the
process. As a result, these methods are becoming increasingly popular in modern theoretical cell biology1,24,40 .

a) Electronic
b) Electronic

mail: abernats@uvm.edu
mail: eguseva@tutanota.com

Additionally, they have more physical rigor compared
to the empirical law of mass action ordinary differential
equations18 .
It is worth mentioning that replacing molecules with
any other interacting agents does not invalidate the approach, as long as certain conditions are met. This makes
it suitable for modeling many non-chemical systems in areas such as population ecology29,37 , evolution theory13,26 ,
immunology11 , epidemiology5,25,27 , sociology12,39 , game
theory30 , economics23,31 , robotics3,4 and information
technology22 .
The majority of contemporary SSAs are based on
the algorithm by Gillespie16–18 known as direct method
(DM). It describes a system with a finite state space undergoing a continuous-time Markov process. The time
spent in every state is distributed exponentially; future
behavior of the system depends only on its current state.
Gillespie has shown18 that for chemical systems these assumptions hold if the system is well-stirred and molecular
velocities follow Maxwell-Boltzmann distribution. Because of DM’s good physical basis its solutions are as
accurate or more accurate than the law of mass equations, which isn’t necessarily a good predictor of mean
values of molecular populations17 .
Time complexity of DM is linear in the number of reactions. For highly coupled systems that means that the
run-time grows as a square of the number of species.
This poses a problem for many models in systems biology which describe systems with a multitude of molecular species connected by complex interaction networks.
Storage complexity of the DM is linear in the number of
possible reactions.
Many alternative SSAs which improve the time complexity of DM (see Table I for an incomplete list) were
developed. Approximate SSAs make such improvements
at the cost of introducing additional approximately sat-

2
isfied assumptions, while exact SSAs achieve better runtime by performing a faster computation that is equivalent to the one performed by DM. Exact methods have
been developed with time per reaction that is linear34 or,
for sparsely connected networks, even constant35 in the
number of molecular species.
One deficiency shared by most SSAs is that they operate with a static list of all possible reactions and species,
while for some systems of interest maintaining such a list
is not possible. An example of such system is heterogeneous polymers undergoing random polymerization, an
object of interest for the researchers of prebiotic polymerization and early evolutionary processes. The number of
heteropolymers that can be produced by such a process
is infinite. Even if we only consider species of length up
to L, we have to deal with O(pL ) species, where p is the
number of monomer species. This causes both time and
storage complexity of stochastic models to grow exponentially with L, which limits the studies to very short
polymers.
However, if an algorithm can maintain a dynamic list
of molecular species and interactions, the complexity of
modeling such systems can be substantially reduced. In
our example, the set of possible species is so large that
even for moderate cutoff L the vast majority of all possible species will have a population of zero. If at any time
step a specie has a population of zero, no reaction can
happen in which this specie is a reagent. Therefore, the
reactions involving such a specie need not be tracked until some reaction occurs in which the specie is a product.
In addition, many studies are concerned with emergent phenomena, which often involve the system exhibiting some dynamical pattern, e.g. sitting in a stable or
chaotic attractor. Such phenomena tend to only involve a
subset of possible states and transitions and not explore
the space of all possibilities uniformly. This may limit
the species and reactions involved in the phenomenon to
a small subset of all possible species and reactions. For
example, in the study by Guseva et al.20 effective populations of all sequences are ∝ 103 , while the number of all
the possible species is ∝ 107 . The model is tightly coupled, thus time costs are 104 times higher and memory
costs are 108 times higher than they could be if only the
final subset of reactions and species was considered and
the SSA retained the time and storage complexity of the
best SSAs for static lists (e.g. PDM34 ).
However, not only this subset is unknown prior to the
experiment, but the transient to the behavior of interest may involve many more species and reactions than
the behavior itself. Thus, the simulation cannot be constrained to use a smaller list of species and reactions if
such list is static.
One class of systems in which reactions cannot be
listed occurs in solid state chemistry. To simulate those,
Henkelman and Jonsson21 developed an algorithm in
which the list of possible reactions is generated on the
fly and used in the standard Gillespie’s algorithm. This
approach was formulated as a part of simulation frame-

work specific to solid state chemistry. It remained obscure outside of the solid state chemistry community and
was rediscovered independently by authors of the present
work early in the course of its preparation.
Here we present extendable partial-propensity direct
method (EPDM): a general purpose, exact SSA in which
species and reactions can be added and removed on the
fly. Unlike Henkelman and Jonsson’s approach, our algorithm is based on the partial-propensity direct method
(PDM) by Ramaswamy et al.34 and retains its linear time
complexity in the number of species. Storage complexity
is linear requirements in the number of reactions.
Only the species with nonzero populations and the reactions involving them are recorded and factored into
the complexity. This reduces time complexity of executing one reaction by a factor of Ntot /N where Ntot is the
number of all possible species and N is the number of
species with nonzero population at the time of the reaction. Storage complexity is reduced by the square of that
factor for densely connected reaction networks.
We test the performance of a C++ implementation of
our algorithm1 for two chemical systems, investigate scaling and compare the performance with two other SSAs.
The results confirm our predictions regarding the algorithm’s complexity.
• Direct method (DM)16,17
• First reaction method (FRM) 16
• Gibson-Bruick’s next-reaction method
(NRM)15
• Optimized direct method (ODM)10
• Sorting direct method (SDM)28
• Partial propensity direct method (PDM)34
Approximate • τ -leaping7,9,19,32
• kα -leaping19
• Implicit τ -leaping36
• The slow-scale method8
• R-leaping2
• L-leap33
• K-leap6
Exact

TABLE I. List of some exact and approximate stochastic simulation algorithms

II.

DIRECT STOCHASTIC ALGORITHMS NOW

A.

Gillespie Algorithm

Being an exact SSA, our algorithm performs an optimized version of the same basic computation as the
Gillespie’s DM14,17 . DM is based on the following observation:

1

https://github.com/abernatskiy/epdm

3
Probability that any particular interaction of
agents occurs within a small period of time
is determined by and proportional to a product of a rate constant specific to this interaction and a number of distinct combinations
of agents required for the interaction.

At every step the algorithm looks through the list of
all M possible reactions. Therefore, the time it takes
to process one reaction (i.e. perform Steps 3 and 4) is
proportional to M :

Since we developed EPDM with chemical applications in
minds, we will hereafter refer to agents of all kinds as
molecules. However, as long as the observation above
holds, the method is valid for any other kind of agent
(see Introduction).
The algorithm starts with initialization (Step 1) of
the types and numbers of all the molecules initially
present in the system, reaction rates and the random
number generator. Then propensities of all reactions are
computed (Step 2). Propensity aµ of a reaction Rµ
(µ ∈ {1, 2, ..., M }) is proportional to its reaction rate
cµ :

There must be a record for every reaction, so the space
complexity is also O(M ).

aµ = hµ cµ .

treac = O(M ).

(1)

Here, hµ is the number of distinct molecular reactant
combinations for reaction Rµ at the current time step,
a combinatorial function which depends on the reaction
type and the numbers of molecules of all reactant types16 .
Total propensity is the sum of propensities of all reactions:
a=

M
X

aµ .

(2)

µ=1

The next step (Step 3), called Monte Carlo step
or sampling step, is the source of stochasticity. Two
real-valued, uniformly distributed random numbers from
[0, 1) are generated. The first (r1 ) is used to compute
time to next reaction τ :
 
1
1
τ = ln
.
(3)
a
r1
The second one (r2 ) determines which reaction occurs
during the next time step τ . The j-th reaction occurs if
j−1
X
µ=1

aj 6 ar2 <

j
X

B.

Partial propensity methods

For many systems it is valid to neglect reactions which
involve more than two molecules or agents. This premise
allows for a class of partial-propensity direct methods.
The method described in the present work is among
those.
If α is the maximum number of reactions which may
happen between any pair of the reagent species then M =
O(αN 2 ). Then the expression for the time complexity of
DM (5) becomes quadratic in N :
treac = O(αN 2 ).

(4)

µ=1

The next step is update (Step 4): simulation time is
increased by τ generated at Step 3, molecules counts are
updated using the stoichiometric numbers of the sampled
reaction and propensities are updated in accordance with
the new molecular counts.
The last step is iteration: go back to Step 3 unless some
termination condition is met. Termination should occur
if no further reactions are possible (i.e. when the total
propensity a = 0). Optional termination conditions may
include reaching a certain simulation time, performing
a given number of reactions, reaching some steady state
etc.

(6)

Partial-propensity direct method (PDM) and sorting
partial-propensity direct method (SPDM)34 improve this
bound to O(αN ) by associating each reaction with one of
the involved reagents and sampling the reactions in two
stages. In the first stage the first reactant specie of the
reaction to occur is determined; this takes O(N ) operations. The second stage determines the second specie
and a particular reaction to occur; this takes O(αN ) operations, and the total complexity of the sampling step
adds up to O(αN ).
Since any specie can be involved in at most O(αN ) unimolecular or bi-molecular reactions, only O(αN ) values
have to be updated when the molecular counts change.
This enables PDM and SPDM to perform the update step
without worsening the time complexity of the sampling
step.
The final time complexity of these algorithms
treac = O(αN )

aj .

(5)

(7)

holds irrespective of the degree of coupling of the reaction
network. The number of records required for sampling is
still O(M ) = O(αN 2 ).
Note that for sparsely coupled reaction networks, time
complexity can be improved further to O(1)35 ; here we
won’t concern ourselves with such reaction networks.

III.
A.

ALGORITHM DESCRIPTION
Reaction grouping

Similar to PDM and SPDM34 , our method relies on
splitting the reactions into groups associated with one

4
of the reagents. To simplify this procedure, we reformulate all reactions with up to two reagents as bimolecular
by introducing the virtual void specie ∅ which always
has a molecular count of 1. It can interact with other,
real species and itself, however its stoichiometry is always zero. With these new assumptions, unimolecular
reactions are reformulated as follows:
Si → products ⇒
1 Si + 0 ∅ → products.

(8)

Source reactions are also reformulated:
∅ → products ⇒
0 ∅ + 0 ∅ → products.

(9)

Our algorithm keeps a list of species existing in the
system to which entries can be added. When a new
specie Si is added to the list, the algorithm considers
every specie Sj in the updated list. For each (unordered)
pairing {Si , Sj } a list of possible reactions is generated.
If the list is not empty, it is associated with Sj , the specie
that has been added to the list earlier unless it coincides
with Si . If Si = Sj , the list is associated with Si .
For example, the first step of the initialization stage
involves adding the first element to the list of species
which is always the void specie ∅. At this point there
are no species in the list aside from ∅, so the algorithm
checks which kinds of reactions may happen between ∅
and itself. Due to the reformulation (9) this will involve
all source reactions. Their list will be generated and associated with the newly added specie ∅.
Another example. Suppose the list of known species
is [∅, S1 ] and we’re adding a specie S2 which reacts with
S1 , itself and also participates in some unimolecular reactions. The list becomes [∅, S1 , S2 ] and the algorithm
proceeds to pair up S2 with every specie in it and generate reaction lists.
Due to the reformulation (8) it will find all the unimolecular reactions for the pair {S2 , ∅}. The list of these
reaction will be associated with the previously known
specie of the pair, ∅.
When considering the pair {S2 , S1 } it will find all the
reactions between S1 and S2 and associate them with S1 .
Finally, it will consider a pair {S2 , S2 } and find its
reactions with itself. Since there are no previously known
species in the pair, the reactions will be associated with
S2 .

B.

Propensities generation

molecules full propensities aµ are
aµ = ni nj cµ for a bimolecular reaction
of distinct species i 6= j, Si + Sj → products,
1
aµ = ni (ni − 1)cµ for a bimolecular reaction
2
(10)
of identical species 2Si → products,
aµ = ni cµ for a unimolecular reaction
Si → products,
aµ = cµ for a source reaction ∅ → products.
Since the void specie ∅ has a fixed population of 1 we can
use the formula for the unimolecular reactions for source
reactions as well: aµ = 1 · cµ for ∅ → products. Then for
all reactions we can define partial propensity w.r.t. the
reactant specie Si as
πµ(i) ≡ aµ /ni .

For the reactions with up to two reagents, partial propensities are34
πµ(i) = nj cµ for a bimolecular reaction
of distinct species i 6= j, Si + Sj → products,
1
πµ(i) = (ni − 1)cµ for a bimolecular reaction
2
(12)
of identical species 2Si → products,
πµ(i) = cµ for a unimolecular reaction
Si → products,
πµ(∅)

= cµ for a source reaction ∅ → products.

Suppose some specie Sj is added to the list of known
species as described in section III A. For every specie Si
in the updated list we generate a list of possible reactions
between Si and Sj and associate it with Si . For every
reaction Rijk in the list, we compute its rate constant cijk
(i)
and its partial propensity πijk w.r.t. Si using formulas
(12).
We also keep some sums of propensities to facilitate
the sampling. Given a list of reactions between Si and
(i)
Sj associated with Si , we define Ψij as
(i)
Ψij

=

αij
X

(i)

πijk ,

(13)

k=1

where αij is the number of reactions possible between Si
and Sj .
(i)
Λi is the sum of partial propensities of all reactions
associated with Si :
(i)

Λi =
As the reactions are generated and associated with
species we also compute and store some propensities.
For all reactions Rµ with up to two participating

(11)

mi
X

(i)

Ψij .

(14)

j=1

Here, mi is the number of lists of reactions with other
species associated with Si .

5
Σi is the full propensity of all reactions associated with
Si :
(i)

Σi = ni Λi .

(15)

and a is defined as a total full propensity of all reactions
in the system:
a=

N
X

Σi .

(16)

i=1

Relations, with as many elements as there are species
with which the first specie has any reactions. On the
other hand, the Population added last will not be an
owner of any Relations and its list of those will be empty.
In this case, all information about the specie’s reactions
is owned by other species and only available within the
Population through the list of RelationAddresses.
A Relation structure owned by Si holds a linked list of
all possible Reactions between Si and Sj and their total
(i)
partial propensity Ψij (see eq. (13)). Reactions store
(i)

C.

Data model

Similar to Steps 3 and 4 in Gillespie’s DM (see section
II A), execution of each reaction in our algorithm involves
two steps: sampling and updating. We’ll refer to the data
used in the sampling process as primary and call all the
data which is used only for updating auxiliary.
To minimize the overhead associated with adding and
removing the data we use a hierarchical linked list (multilevel linked list of linked lists). We define several data
structure types (see table II), many of which include
lists of instances of other types. Complete data model
is shown in Figure 1.
Top level list is stored within a structure of type
TotalPopulation. Aside from the list, the structure
holds the data describing the system as a whole: total
propensity a (see (16)), current simulation time t and
the random generator state.
Elements of the top level list are structures of type
Population, each of which holds the information related
to the molecular population of a particular specie Si . The
information about the specie itself, to which we will in
the context of its Population refer as the owner specie,
is represented by a member object of user-defined class
Specie (see note 2 at the end of this section), containing
a unique string specie identifier IDi as a member variable. Population also contains the number of molecules
in the population ni , total propensity of the population
(i)
Σi and total partial propensity Λi w.r.t. Si of all reactions associated with Si (see eqs. (15) and (14)).
The Population structures are appended to the top
level list sequentially over the course of the algorithm
execution (see section III A). Each of them holds, in addition to the data mentioned above, two linked lists: of
structures of type Relation and of structures of type
RelationAddress. The former stores the list of all reactions possible between the owner specie and some other
specie that has been added later in the course of the algorithm’s operation. RelationAddress structures hold the
references used to access Relations in which the owner
specie participates, but not as an owner. Those include
all the species added before the owner.
It can be observed that the Population which has
been added first will necessarily has an empty list
of RelationAddresses and a potentially large list of

the reaction rate cijk , partial propensity πijk , a table of
stoichiometric coefficients and IDs of all participating
species, including products. Additionally, each of them
stores an auxiliary reference to the RelationAddress
structure pointing to this Relation.
RelationAddress is an auxiliary structure holding references to a Relation in which the owner specie participates without owning it, and to the Population of the
specie which does own the Relation. It also contains the
references to itself and to the list holding it.
Notes:
1. We will say that a TotalPopulation is valid iff
all variables mentioned in eqs. (12–16) satisfy
these equations and all references mentioned in the
present section are valid and point where intended.
2. Unique specie string ID format and Specie objects
are user-defined and must be freely convertible between each other. The utility of having a separate
Specie class lies in it having a user-defined method
reactions() : (self, Specie) → ListOf Reactions
which produces a list of all reactions possible between the specie described by a caller object and a
specie described by the the argument object. For
some systems, this computation can be made much
faster if some auxiliary data can be kept in the
structure describing a specie. String IDs on the
other hand are intended as memory efficient representations of specie data which are used when this
computation may not be needed, e.g. to represent
reaction product not yet present in the system or
for specie comparison.
Our method does not rely on this separation other
than as a means of optimization.
3. Additionally, in our implementation we provide
a global structure to hold the parameters of the
system which may influence the behavior of the
reactions() method of class Specie.
D.

Adding a Population

Suppose that we have a valid TotalPopulation with
N − 1 Populations. To add a Population based
on a specie ID and a molecular count n, we follow
the method described in section III A (see also algorithm 1). We begin by appending a Population pN to

6
Data type

Containers
TotalPopulation Linked
list
Populations

Population

Relation

Linked
list
of
Relations,
linked
list
of
RelationAddresses
Linked
list
of
Reactions

RelationAddress –

Reaction

Specie

of

Necessary members
Scalars
References
Total propensity –
a (16), current
time t, random
generator state
Specie S, molec- –
ular count n,
propensity sums
Σ (15) and Λ (14)
Partial propensity To
sum Ψ (13)
RelationAddress
pointing to this
Relation
–
To a Relation,
its
owner’s
Population, list
containing
this
RelationAddress
and itself
Rate constant c, –
partial propensity
π w.r.t. the owner
specie (12)

Array
of
pairs
(ID, σ) of IDs
and stoichiometries
of
participating
species
–
Unique, compact –
specie identifier
ID

Significance
Functions
–
Represents the whole system being modeled

–

Represents a population
molecules of a specie S

–

Data on all reactions possible
between a pair of species, to
be stored withing the list of an
owner specie
Reference to a Relation to be
kept by a non-owner specie, useful in propensity updates and
population deletions

–

–

of

Represents a reaction

Specie:: An extended specie represenreactions tation capable of keeping ex(Specie) tra information to generate
lists of possible reactions with
reactions() method quickly

TABLE II. List of data structure types used by the algorithm

TotalPopulation’s list. pN ’s Specie structure SN is
converted from ID, its molecular count nN ← n and
(N )
propensities ΛN ← 0, ΣN ← 0.
Next, for each Population pi in TotalPopulation’s
updated list (including the newly added pN ) we compute
a list of all possible reactions between SN and pi ’s Specie
Si . If the list is not empty, it is then converted into
Relation ρiN and stored at pi ’s list of those. To do
(i)
the conversion, we compute partial propensities πiN k and
(i)
ΨiN using eqs. (12) and (13). We also store copies of Si ’s
ID at every Reaction in the list and at the Relation ρiN
itself, to keep track of the specie w.r.t. which the current
partial propensities are computed. ρiN ’s reference to the
RelationAddress structure is invalid at this point.
After appending the Relation we update pi ’s propen(i)
(i)
(i)
(i)
sities: Λi ← Λi + ΨiN , Σi ← ni Λi .
We then proceed to create a RelationAddress structure RAN i pointing to ρiN . A blank structure is appended to pN ’s list. We take references to ρiN , pi , RAN i
and pN ’s list of RelationAddresses and save them at
RAN i .
Then, we take the reference to RAN i and store it at
ρiN . At this point, all references in the whole structure
are valid and correct.
Processing each pi in the TotalPopulation’s list takes
O(α), so the whole procedure up to this point takes

O(αN ). At this point we need to update the total
propensity of the system a by recomputing it, which takes
O(N )2 .
The total number of operations it takes to add one
Population is O(αN ). The operation preserves the validity of the data structure.
E.

Initialization stage

To initialize the data structure, we make a
TotalPopulation with an empty list of Populations,
initialize a ← 0, t ← 0 and the random generator with
a seed. The only variable that needs to take a particular value for the TotalPopulation to be valid is
a and it has the correct value of 0; therefore, it is a
valid TotalPopulation. We build the data structure by
adding the initial populations to this structure as described in section III D. The resulting structure is valid
because we only used validity-preserving operatins.
Adding every population takes O(αN ) operations and
it must be repeated N times. This brings the complexity
of the initialization step to O(αN 2 ).

2

It can be done during the iteration in O(1), but we chose to
recompute it for improved numerical accuracy.

7
TotalPopulation tp

a, t, RNG
Population pi

Population pj

ni, Σi, Λi

nj, Σj, Λj

Specie Si

Specie Sj

IDi

IDj

Relation ρij

Relation ρjl

Ψjl

Ψij

From Relation ρoi with
some Specie So ≠ Sj

σj(Rjlk') IDj
σl(Rjlk') IDl
σp1(Rjlk') IDp1
...
σpm(Rjlk')IDpm

products

RelationAddress RAio

To Relation
Address RAlj

Reaction Rjlk'
cjlk', πjlk'(j)

reactants

products

σi(Rijk) IDi
σj(Rijk) IDj
σp1(Rijk) IDp1
...
σpm(Rijk)IDpm

reactants

Reaction Rijk
cijk, πijk(i)

From Relation
Address RAlj held
by the Population
of Sl

RelationAddress RAji

To Population of
the Specie So
To Relation of
the Specie So to Si

FIG. 1. EPDM data model. Circles with arrows denote references. References shown in black form the hierarchical linked
list; in implementation, many of those are hidden within the std::list container. Auxiliary references added for dynamic
updates are shown in blue. Letters and boxes within the boxes indicate non-reference member variables. Crossed empty circle
represents a null reference.

F.

Deleting a Population

Our algorithm is designed to keep track only of the
species with a nonzero molecular count and reactions involving them. To accomplish that, we use the addition
operations described in the previous section and deletion
operation described here (see also algorithm 2). Deletion
is only ever applied to Populations of species with zero
molecules. All propensities of reactions involving such
species are zeros; this enables us to simplify the procedure.
To remove a Population pi , we begin by removing
all RelationAddress structures pointing at Relations
in pi ’s list. Each Relation ρij in pi ’s list contains a
reference to the RelationAddress structure pointing at

it, RAji , which in turn contains a reference to itself and
to a list holding it. We use those to remove each RAji
from its list. Since Si can be involved in at most N
relations, this step takes O(N ) operations.
Next, we remove all Relations in which Si participates, but which are owned by other species. For each
RelationAddress RAij we follow its references to the
Population pj of the other specie and to its Relation
ρji with Si . We use those to remove ρji from pj ’s list.
This step also takes O(N ) operations.
Finally, we delete the whole Population structure pi
from TotalPopulation’s list. We recursively remove
all the structures in its lists, which takes O(N ) operations for the list of RelationAddresses and O(αN )
operations for the list of Relations. The resulting

8
G.

Sampling stage

Algorithm 1 Adding a Population
1: function AddPopulation(tp, ID, n)
2:
Append a new Population pN to tp’s list
3:
Convert ID to a Specie object S
4:
SN ← S, nN ← n
(N )
5:
ΛN ← 0, ΣN ← 0
6:
for all Populations pi in tp’s list do
7:
ListOf Reactions ← SN .reactions(Si )
8:
if ListOf Reactions is empty then
9:
continue
10:
end if
11:
Append a new Relation ρiN to pi ’s list
12:
Append all Reactions in ListOf Reactions
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:

to ρiN ’s list
(i)
ΨiN ← 0
for all Reactions RiN k in ρiN ’s list do
(i)
Compute πiN k using (12)
(i)
(i)
(i)
ΨiN ← ΨiN + πiN k
end for
(i)
(i)
(i)
Λi ← Λi + ΨiN
Append a new RelationAddress structure RAN i
to pN ’s list
Store references to ρiN , pi , RAN i and pN ’s list of
RelationAddresses at RAN i
Store reference to RAN i at ρiN
end for
Recompute a using (16)
end function

Similarly to DM16 and PDM34 , our algorithm simulates the system by randomly sampling time to the next
reaction and the reaction itself with certain distributions.
Here we describe how it happens in our algorithm.
We begin by generating two random numbers r1 and
r2 . The first one is used to compute the time to the next
reaction τ exactly as in DM and PDM (see eq. (3)). The
second random number r2 is used to sample the reaction
similarly to how its done in PDM34 .
The sampling process has three stages (see algorithm
3). During the first stage (lines 2-6) we determine the
first specie participating in the reaction to happen. To
this end, we go through the list of Populations p1 ...pN
until the following condition is satisfied:
J
X

J
X

Σ i + nJ

K
X

(J)

ΨJi 6 ar2 <

4:
5:
6:
7:
8:
9:
10:
11:
12:

pointing at it, RAji
Use the reference at RAji to itself and to the
list holding it to remove it from the list
end for
for all RelationAddress RAij in pi ’s list do
Follow RAij ’s reference to the Relation ρji
and to the Population pj owning it
Remove ρji from pj ’s list
end for
Delete all the date in the pi structure
Remove pi from tp’s list
end function

(17)

J
X

Σi + n J

K+1
X

i=1

i=1

(J)

ΨJi (18)

i=1

During the third stage (lines 12-16), we determine
which of the reactions possible between the two species
is going to happen. We go through ρJK ’s list of the
Reactions, looking for a reaction RJKL such that
J
X

Σi + n J

i=1

K
X

(J)

ΨJi + nJ

i=1

<

J
X

Σ i + nJ

i=1

Algorithm 2 Deleting a Population

Σi .

i=1

The second stage (lines 7-11) involves finding the second reactant. We look for a Relation ρJK among those
attached to the Population pJ for which the following
condition holds:

The final complexity of the deletion operation is
O(αN ).

1: function DeletePopulation(tp, pi )
2:
for all Relations ρij in pi ’s list do
3:
Follow ρij ’s reference to RelationAddress

J+1
X

i=1

i=1

TotalPopulation is valid since all the propensities of
the reacitons involving Si are zeros and the remaining
propensity sums has not changed; it also contains no invalid references.

Σi 6 ar2 <

L
X

(J)

πJKi 6 ar2 <

i=1
K+1
X
i=1

(J)
ΨJi

+ nJ

L+1
X

(19)
(J)
πJKi .

i=1

Note how equations (17) and (18–19) are similar, but
their implementation in the pseudocode (algorithm 3) is
different. This design minimizes sampling errors due to
floating point representation, making the first stage exact.
The resulting reaction sampling finds a reaction in exactly the same manner as equation (4) does. However,
the first and the second stages of this sampling process
take O(N ) steps and the third step takes O(α) steps, resulting in a total time complexity of O(N + α). Using
(4) directly requires O(M ) steps16 , which is O(αN 2 ) for
densely connected reaction networks.

H.

Updating stage

When the reaction to occur RJKL is known, our algorithm proceeds to update the data to reflect the changes

9

Algorithm 3 Reaction sampling

Algorithm 4 Post-sampling data update

1: function SampleReaction(tp, r2 )
2:
s1 ← 0, s2 ← 0
3:
while ar2 > s1 do
4:
Get a new Population pi from tp’s list
5:
s2 ← s1 , s1 ← s1 + Σi
6:
end while
7:
g ← (ar2 − s2 )/ni
8:
while g > 0 do
9:
Get a new Relation ρij from pi ’s list
10:
11:
12:
13:
14:
15:
16:
17:
18:

1: function UpdateExistingPopulations(tp, RIJK )
2:
for all species Si participating in RIJK do
3:
Search for a Population of specie Si in tp’s list
4:
5:
6:
7:

(i)

g ← g − Ψij
end while
(i)
g ← g + Ψij
while g > 0 do
Get a new Reaction Rijk from ρij ’s list
(i)
g ← g − πijk
end while
return Rijk
end function

8:

←

(j)
Λj

−

(j)

(j)
Ψji ,

11:
12:
13:
14:

Ψji ← Ψji + πjik
end for
(j)
(j)
(j)
Λj ← Λj + Ψji

(j)

(j)

(j)

(j)

15:
Σj ← nj · Λj
16:
end for
17:
else
18:
AddPopulation(tp, IDi , σi (RIJK ))
19:
end if
20:
end for
21:
for all Populations pi in tp’s list do
22:
if ni == 0 then
23:
DeletePopulation(tp, pi )
24:
end if
25:
end for
26: end function

(20)

Stoichiometric coefficients are negative for reagents, so
their molecular counts may become zero after this step.
After updating the molecular count, we also update all
partial and total propensities which depend on it. From
every RelationAddress RAij in pi ’s list we obtain references to a Relation ρji in which Si participates and
to the Population pj of its owner specie. For ρji we
(j)
(j)
recompute all πjik and Ψji from scratch using formulas
(12) and (13). For pj , we update the propensity sums as
follows:
(j)
Λj

(j)

Ψji ← 0
for all Reactions Rjik in ρji ’s list do
(j)
Recompute πjik using (12)

9:
10:

in species’ populations and propensities (see also algorithm 4). For every specie Si involved in RJKL we read
its IDi from the array stored at RJKL . We run a sequential search for this specie’s Population pi over the
list kept in TotalPopulation. If the specie’s Population
is found, its molecular count ni is updated using the stoichiomethic coefficient σi (RJKL ):
ni ← ni + σi (RJKL ).

using its IDi
if Si has a Population pi in tp’s list then
ni ← ni + σi (RIJK )
for all RelationAddresses RAij in pi ’s list do
Follow the references at RAij to get
the Relation ρji and its owner
Population pj
(j)
(j)
(j)
Λj ← Λj − Ψji

(21)

Σj ← nj Λj .
Since each of the species involved in RJKL may be
involved in O(αN ) reaction, updating the structure in
this way takes a total of O(αN ) operations.
If the specie Si is a product, its Population may not
exist yet. In this case we add a new Population of Si
using its IDi as described in section III D. The molecular count of the newly added specie is its stoichiometric
coefficient in RJKL , σi (RJKL ). Additions take O(αN )
operations.
When we’re done updating the existing Populations
and adding the new ones, we iterate through the list of
Populations again and delete the ones with a molecular

count of zero as described in section III F. The deletion
takes O(αN ).
Finally, we recompute the total propensity of the system a using (16).

I.

Summary of EPDM

EPDM is an SSA which only maintains the data about
the species with nonzero molecular count (see algorithm
5). This ensures that the number of tracked molecular
species N is as low as possible.
Our algorithm uses a data structure described in section III C. The structure holds one entry for each possible
reaction, bringing storage requirements of our algorithm
to O(αN 2 ).
The data structure is initialized by constructing it to
be empty, then adding the specie data for every specie
initially present in the system. The process is described
in section III E and takes O(αN 2 ) operations.
Each step of the simulation executes a single reaction.
It is composed of a sampling step (see section III G) and
an updating step (section III H). Each of these takes αN
operations, so the total number of operations needed to
simulate one reaction is also O(αN ).

10
IV.

BENCHMARKS

Algorithm 5 EPDM overview
1: . Initialization stage, see section III E
2: Create TotalPopulation tp based on initial conditions
3: while a > 0 and termination conditions not met do
4:
. Sampling stage, see section III G
5:
Generate two random numbers r1 and r2
6:
Compute τ using (3)
7:
R ←SampleReaction(tp, r2 )
8:
. Updating stage, see section III H
9:
for all species Si involved in R do
10:
Search tp’s list for Population of Si
11:
if a Population pi of Si has been found then
12:
Update Si ’s molecular count ni
13:
Update all partial propensities which

We benchmark performance of EPDM against two direct methods: PDM343 and DM384
Because our model was designed with complex systems
in mind we studied performance only for strongly coupled
systems. In both systems every specie can interact with
every other specie in a unique way, ensuring that the
number of reactions is
M=

When some reaction produces any number of molecules
of a previously unknown specie, specie data is added as
described in section III D. To generate the reactions dynamically, a user-defined function reactions() is used
which takes two species and produces a list of reactions
between them, complete with rates.
When any specie’s molecular count reaches zero, its
data is pruned from the structure as described in section
III F.

J.

Implementation

We implemented our algorithm as a C++ framework.
The user must define a class Specie with a constructor
from a string ID which must save the string into the
member variable m id. The class must define a method
Specie::reactions(Specie) returning a list of possible
Reactions between the caller Specie and the argument.
After implementing the class Specie, users can simulate the system. Two stopping criteria are currently
available by default: the algorithm can stop either when
a certain number of reactions have been executed or a
certain simulation time has passed.
A global dictionary with arbitrary parameters loaded
from a configuration file is provided for convenience.
The implementation is currently available for Linux
and Mac OS X. It is tested with GNU gcc 4.9.3 and
GNU make 4.1.
The code is available at https://github.com/
abernatskiy/epdm.

(22)

and

depend on ni
14:
else if no Population of Si has been found then
15:
Add a Population of Si to tp’s list
16:
. details in section III D
17:
end if
18:
end for
19:
Find all Populations with n == 0 in tp’s list
20:
Remove all found Populations
21:
. details in section III F
22: end while

N (N + 1)
2

α = 1.

(23)

Both models are designed in such a way that the total
number of species is preserved throughout the simulation
time.
Models below are designed to measure the performance
of our method. They don’t fully illustrate the power of
the model because they have fixed number molecules,
which is necessary to keep for benchmark. The most
striking performance gain is achieved when listing all the
species isn’t possible in principle or due to computational
costs. For example, in case of realistic polymerization
and autocatalysis model used to study prebiotic polymerization Guseva et al. it was possible to increase the
maximum length of simulated polymers from 12 to 25 by
employing our algorithm. Note that limit of 25 wasn’t
due to restrictions of our algorithm, but due to necessity to calculate minimum energy folding configuration
of every chain, which is an NP-complete problem.
A. CPU time of EPDM is linear for a strongly coupled
system

The first model (”colliding particles”) is made to test
how our algorithm performs on chemical systems where
no new molecules are created and no molecules ever disappear. This is a model of a system consisting of colliding
particles of N species. Particles behave like rigid spheres:
they collide, and bounce back without internal changes.
All of the particle species are known in advance, none are
added or removed over the course of the simulation. The
system is defined by the following set of equations:
k

Pi + Pj −
→ Pi + Pj ,

i, j ∈ [1, N ].

(24)

In our simulations we vary number of species N from 10
to 7000. Every specie has a population of 50 molecules.

3
4

We
took
implementation
from
http://mosaic.mpicbg.de/pSSALib/pSSAlib.html
We
took
implementation
from
http://sourceforge.net/projects/stochkit/

11

4

2
1

EPDM
PDM
DM

0
0

2500

N

5000

7500

FIG. 2. CPU time per reaction as a function of the number of
species in the system for ”colliding particles” model (section
IV A) for EPDM, PDM and DM.

Collision rate k is fixed at 0.5 s−1 . Every simulation
runs until 5000 reactions have occurred. For every value
of N , CPU time to simulation completion was measured
10 times and the average time is reported.
Figure 2 shows CPU time it takes per reaction for DM,
PDM and EPDM. DM is clearly quadratic. For any given
N PDM outperforms the EPDM, but both are linear in
time. It is important to note that DM and PDM were
stopped for a relatively low values of N due to excessive
RAM consumption (about ∼ 120GB) by both of the applications, which we suspect was due to implementation
issues (in particular, in XML handling libraries).

B. CPU time stays linear when species are actively
deleted and created

The second test checks if algorithm keeps linear time
when Populations are added and removed from the system. The test system (”particles with color”) consists of
N colliding particles of N types, each of which has an internal property (”color”) that is changed in the collision.
The following equation defines the system:
k

Piα + Pjβ −
→ Piα+1 + Pjβ+1 ,

50

i, j ∈ [1, N ], α, β ∈ [1, Ω]
(25)
Indexes i, j enumerate particle types and run from 1 to
N . Indexes α, β enumerate colors of particles; particle
can have one of Ω colors. During the collision color index
of a participating particle goes up Piα → Piα+1 , until it
reaches maximal index Ω, after which it drops back to 0:
PiΩ → Pi0 .
Since every combination of a particle type and a color
is considered a separate specie, every reaction causes
two species to go extinct and their Populations to be
deleted. It also adds two new species, which requires

t, ms

t, ms

3

adding two new Populations. Thus, the total number
of species simultaneously present in the system is maintained at exactly the same level. Every specie is represented in the system by a single molecule.
To run such a simulation in DM and PDM frameworks
we had to enumerate all N Ω of the possible species. This
slows down the simulation enough to make the comparison impossible beyond a small number of species. The
figure 3 shows how CPU time per reaction depends on
the number of species in EPDM.

25

0
0

2500

N

5000

7500

FIG. 3. CPU time per reaction for EPDM simulation of ”particles with color” (see section IV B).

V.

CONCLUSION AND DISCUSSION

Stochastic simulations are actively used in molecular
and systems biology. The bigger and more complex the
system, the more important the performance of the simulation algorithm becomes. It is also more burdensome
or even impossible to list all the species and reactions for
complex systems. We introduced a general purpose, exact stochastic simulation algorithm which allows to avoid
listing all the possible species and reactions by defining
the general rules governing the system instead. Built
within the partial propensity framework34,35 , our algorithm achieves linear time complexity in the number of
molecular species.
The algorithm has its limitations. First, it is limited to
reactions of maximum of two reactants. In the chemical
system that doesn’t present much of an issue because reactions with three molecules are significantly rarer then
binary reactions and more complex reactions can be represented as sequences of binary reactions. Second, it cannot simulate spatially nonhomogeneous systems. However, as long as all reactions involve no more than two
reactants and the observation from section II A holds, it
is possible to simulate any system for which the set of
reactions between any two species is known.
Benchmarks suggest that our algorithm is slower than
PDM by a constant factor. Therefore, one should use

12
PDM when it is unlikely that a significant proportion of
species will have a molecular count of zero.
Our results suggest that in complex systems (such as
polymerization reactions and networks of intra-cellular
reactions) and in the case of emergent phenomena studies (e.g. origins of life) EPDM can significantly improve the performance of the stochastic simulations. The
software implementation of the algorithm is available as
an open source public repository https://github.com/
abernatskiy/epdm.

VI.

ACKNOWLEDGMENTS

Authors appreciate support from the NSF grant
MCB1344230.
1 Arkin,

A., Ross, J., and McAdams, H. H. (1998). Stochastic
kinetic analysis of developmental pathway bifurcation in phage
lambda-infected Escherichia coli cells. Genetics, 149(4):1633–48.
2 Auger, A., Chatelain, P., and Koumoutsakos, P. (2006). Rleaping: Accelerating the stochastic simulation algorithm by reaction leaps. The Journal of Chemical Physics, 125(8):084103.
3 Berman, S., Halasz, A., Hsieh, M., and Kumar, V. (2009). Optimized Stochastic Policies for Task Allocation in Swarms of
Robots. IEEE Transactions on Robotics, 25(4):927–937.
4 Berman, S., Halasz, A., Kumar, V., and Pratt, S. (2007). BioInspired Group Behaviors for the Deployment of a Swarm of
Robots to Multiple Destinations. In Proceedings 2007 IEEE
International Conference on Robotics and Automation, pages
2318–2323. IEEE.
5 Breban, R., Drake, J. M., Stallknecht, D. E., and Rohani, P.
(2009). The Role of Environmental Transmission in Recurrent Avian Influenza Epidemics. PLoS Computational Biology,
5(4):e1000346.
6 Cai, X. and Xu, Z. (2007). K -leap method for accelerating
stochastic simulation of coupled chemical reactions. Journal of
Chemical Physics, 126(7):1–10.
7 Cao, Y., Gillespie, D. T., and Petzold, L. R. (2005a). Avoiding
negative populations in explicit Poisson tau-leaping. The Journal
of Chemical Physics, 123(5):054104.
8 Cao, Y., Gillespie, D. T., and Petzold, L. R. (2005b). The slowscale stochastic simulation algorithm. The Journal of Chemical
Physics, 122(1):014116.
9 Cao, Y., Gillespie, D. T., and Petzold, L. R. (2006). Efficient
step size selection for the tau-leaping simulation method. The
Journal of Chemical Physics, 124(4):044109.
10 Cao, Y., Li, H., and Petzold, L. (2004). Efficient formulation of
the stochastic simulation algorithm for chemically reacting systems. The Journal of Chemical Physics, 121(9):4059.
11 Caravagna, G., dOnofrio, A., Milazzo, P., and Barbuti, R. (2010).
Tumour suppression by immune system through stochastic oscillations. Journal of Theoretical Biology, 265(3):336–345.
12 Carbone, G. and Giannoccaro, I. (2015). Model of human collective decision-making in complex environments. The European
Physical Journal B, 88(12):339.
13 Dieckmann, U. and Law, R. (1996). The dynamical theory of
coevolution: a derivation from stochastic ecological processes.
Journal of Mathematical Biology, 34(5-6):579–612.
14 Doob, J. L. (1953). Stochastic Processes. John Wiley & Sons,
Inc.;Chapman & Hall, New York, New York, USA.
15 Gibson, M. A. and Bruck, J. (2000). Efficient Exact Stochastic
Simulation of Chemical Systems with Many Species and Many
Channels. The Journal of Physical Chemistry A, 104(9):1876–
1889.

16 Gillespie,

D. T. (1976). A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.
Journal of Computational Physics, 22(4):403–434.
17 Gillespie, D. T. (1977).
Exact stochastic simulation of coupled chemical reactions. The Journal of Physical Chemistry,
81(25):2340–2361.
18 Gillespie, D. T. (1992). A rigorous derivation of the chemical
master equation. Physica A: Statistical Mechanics and its Applications, 188(1-3):404–425.
19 Gillespie, D. T. (2001). Approximate accelerated stochastic simulation of chemically reacting systems. The Journal of Chemical
Physics, 115(4):1716.
20 Guseva, E. A., Zuckermann, R. N., and Dill, K. A. (2016). How
did prebiotic polymers become informational foldamers?
21 Henkelman, G. and Jonsson, H. (2001). Long time scale kinetic Monte Carlo simulations without lattice approximation
and predefined event table. The Journal of Chemical Physics,
115(21):9657.
22 Hillston, J. (2005). Fluid flow approximation of PEPA models. In
Second International Conference on the Quantitative Evaluation
of Systems (QEST’05), pages 33–42. IEEE.
23 Juanico, D. E. (2012). Critical network effect induces business
oscillations in multi-level marketing systems.
24 Karlebach, G. and Shamir, R. (2008). Modelling and analysis of
gene regulatory networks. Nature Reviews Molecular Cell Biology, 9(10):770–780.
25 Keeling, M. J. and Rohani, P. (2008). Modeling infectious diseases in humans and animals. Princeton University Press.
26 Le Galliard, J.-F., Ferrière, R., and Dieckmann, U. (2005).
Adaptive evolution of social traits: origin, trajectories, and correlations of altruism and mobility. The American naturalist,
165(2):206–24.
27 Legrand, J., Grais, R. F., Boelle, P. Y., Valleron, A. J., and
Flahault, A. (2007). Understanding the dynamics of Ebola epidemics. Epidemiology and infection, 135(4):610–21.
28 McCollum, J. M., Peterson, G. D., Cox, C. D., Simpson, M. L.,
and Samatova, N. F. (2006). The sorting direct method for
stochastic simulation of biochemical systems with varying reaction execution behavior. Computational Biology and Chemistry,
30(1):39–49.
29 McKane, A. J. and Newman, T. J. (2005). Predator-prey cycles
from resonant amplification of demographic stochasticity. Physical review letters, 94(21):218102.
30 Mobilia, M. (2012).
Stochastic dynamics of the prisoner’s
dilemma with cooperation facilitators. Physical review. E, Statistical, nonlinear, and soft matter physics, 86(1 Pt 1):011134.
31 Montero, M. (2008). Predator-Prey Model for Stock Market Fluctuations. SSRN Electronic Journal.
32 Peng, X., Zhou, W., and Wang, Y. (2007). Efficient binomial
leap method for simulating chemical kinetics. The Journal of
Chemical Physics, 126(22):224109.
33 Peng, X.-j. and Wang, Y.-f. (2007). L-leap: accelerating the
stochastic simulation of chemically reacting systems. Applied
Mathematics and Mechanics, 28(10):1361–1371.
34 Ramaswamy, R., González-Segredo, N., and Sbalzarini, I. F.
(2009). A new class of highly efficient exact stochastic simulation algorithms for chemical reaction networks. The Journal of
chemical physics, 130(24):244104.
35 Ramaswamy, R. and Sbalzarini, I. F. (2010). A partial-propensity
variant of the composition-rejection stochastic simulation algorithm for chemical reaction networks. Journal of Chemical
Physics, 132(4):1–6.
36 Rathinam, M., Petzold, L. R., Cao, Y., and Gillespie, D. T.
(2003). Stiffness in stochastic chemically reacting systems: The
implicit tau-leaping method. The Journal of Chemical Physics,
119(24):12784.
37 Reichenbach, T., Mobilia, M., and Frey, E. (2007). Mobility promotes and jeopardizes biodiversity in rockpaperscissors games.
Nature, 448(7157):1046–1049.

13
38 Sanft,

K. R., Wu, S., Roh, M., Fu, J., Lim, R. K., and Petzold,
L. R. (2011). StochKit2: software for discrete stochastic simulation of biochemical systems with events. Bioinformatics (Oxford,
England), 27(17):2457–8.
39 Shkarayev, M. S., Schwartz, I. B., and Shaw, L. B. (2013). Recruitment dynamics in adaptive social networks. Journal of

physics. A, Mathematical and theoretical, 46(24):245003.
D. J. (2011). Stochastic Modelling for Systems Biology. CRC Press.

40 Wilkinson,

