A stochastic model of autocatalytic reaction
networks
Alessandro Filisetti1, Roberto Serra2;1, Marco Villani2;1, Rudolf M.
Fuchslin3;1, Norman Packard4;1, Stuart A. Kauman5, and Irene Poli6;1
1 European Centre for Living Technology
Calle del Clero 2940, 30124 Venice, Italy
alessandro.filisetti@ecltech.org
2 Dipartimento di Scienze Sociali, Cognitive e Quantitative
Universita di Modena e Reggio Emilia, via Allegri 9, 42100 Reggio Emilia, Italy
rserra@unimore.it, villani.marco@unimore.it
3 Articial Intelligence Lab Univ. Zurich
Andreasstr. 15, CH-8050 Zurich, Switzerland
Oce: +41 (0)44 635 45 9, Mobile: +41 (0)79 232 74 36
fuchslin@ifi.uzh.ch
4 ProtoLife Inc
57 Post St. Suite 513, San Francisco, CA 94104
n@protolife.com
5 Departments of Biochemistry and Mathematics
University of Vermont, Burlington, VT 05405
Stuart.Kauffman@uvm.edu
6 Dipartimento di Statistica, Universita Ca' Foscari,
San Giobbe - Cannaregio 873, 30121 Venezia, Italy
irenpoli@unive.it
Abstract. Autocatalytic networks are widespread in nature, but they
are dicult to create or to reproduce in laboratory. There are however
several models of coupled reactions which describe a phase transition to
an autocatalytic cycle when a certain level of heterogeneity in the com-
position of the chemical soup is reached, so it is interesting to understand
why these phenomena are not easily achieved in the laboratory. For this
purpose we introduce here a model, inspired by a previous one by Kau-
man, tailored for the study of such properties. In particular, we take into
account the stochastic nature of the dynamics of interacting molecules,
in the case of a well stirred tank reactor. We describe the model and
we analyse its behaviour under dierent circumstances. In particular,
the onset of an autocatalytic set is studied as the feed is varied, and its
stability is analysed.
1 Introduction
There are many examples of autocatalytic networks in present day biological
systems, which are the outcome of billion of years of evolution. However, if one
is interested either in the problem of the origin of life or in the design of articial

2 Filisetti et al.
protocells [3, 7, 16{18], it is important to study the generic properties of these
systems.
We introduce here a model tailored for the study of such properties. Some mod-
els of this kind have been developed in the past, including the well-known works
by Dyson [4], Eigen and Schuster [5], Kauman [12], Jain and Khrishna [10] and
Kaneko [11].
The major hypotheses about the physical nature of the collectively replicating
molecules and of their interactions can be grouped in two classes: one which is
typical of gene-rst scenarios and is inspired by nucleic acid template matching,
where a linear polymer drives the synthesis of its complementary strand, and one
which is inspired by protein-rst scenarios, where one supposes that a molecule
can catalyse reactions involving other molecules without imposing the constraint
of complementary matching.
Note that nowadays there is a wide consensus that the appearance of collectively
self-replicating sets of molecules is not per se \life", whose properties require
also the existence of a container which separates the \living system" from its
environment, a coupling between the dynamics of the replicating molecules and
the growth and division of the container [17, 7, 3]. Once such coupling has been
achieved, it has been proven that the rates of the two processes (duplication of
the genetic molecules and of the container) tend to spontaneously synchronise
through successive divisions [15]. In any case, the appearance of a set of collec-
tively self-replicating molecules is a crucial step both for the origin of life and
for that of an eective protocell.
Note that, in spite of their important dierences, all the above mentioned mod-
els describe a phase transition in the system behaviour, so that when a certain
condition is met an autocatalytic set appears. On the other hand it is observed
in experiments that obtaining a set of molecules endowed with such a property
is (at least) a dicult task, and this seems of course to contradict the apparent
universality of the theoretical behaviour.
There are essentially three ways out of this impasse. One (which is preferred by
some experimentalists) claims that the simplifying assumptions of these models
make them too unrealistic to be able to say interesting things about the real
world and its molecules. The second one (advocated by some theorists) is that
the true experiment has never been done, because the threshold for collective
autocatalysis has never been reached in experiments with random molecules.
We take here a third position, which is that of trying to analyse which phenom-
ena can aect the behaviour of the models, and to see whether the transition
is retained when they are taken into account and, in case, to see whether the
threshold is severely changed.
In particular, in this paper we wish to take rst into account the stochastic na-
ture of the dynamics of interacting molecules, which may lead to counterintuitive
behaviours if the system is close to an instability point and if the numerosity of
some of its species is low.
Our starting point will be one of the above-mentioned models, i.e. that proposed
by Kauman [12], which is summarised in section 2. Brie
y, the molecules are

A stochastic model of autocatalytic reaction networks 3
linear chains of monomers taken from a nite alphabet. There are two possible
reactions, namely condensation (two polymers are joined forming a longer one)
and cleavage (a polymer gives rise to two by splitting at a certain point). It is
assumed that these reactions occur at a negligible rate unless they are catalysed,
and it is also assumed that any molecule has a small but nite chance to catalyse
a given reaction. By enumerating all the possible reactions and molecules, Kau-
man came to the conclusion that, provided that there are enough molecules in
the initial set, a connected component will appear in the reaction graph, this cor-
responding to the presence of (at least) an autocatalytic set (ACS from now on).
In section 2 we show that, while the probability of assembling a single specic
polymer is a rapidly decreasing function of its length, the number of possible
polymers increases very fast with chain length, so that the formation of long
polymers is favoured.
The reasoning of Kauman about the transition is compelling, but the simplica-
tions were drastic. In particular, in his original work Kauman did not consider
selection nor the concentration of the molecules, but focused his attention on
the graph of the reactions which are possible among molecules which can be
synthesised. Indeed, no extinction, no small-number eect.
A further step in this direction was taken by Farmer and others [6, 1], who
switched to a representation in terms of dierential equations, thereby taking
the concentrations into account. Moreover, they modelled a chemostat, thereby
providing a kind of selective pressure. Their results are interesting, and will be
brie
y summarised in section 2. However, the use of a continuous formalism
does not allow to take stochastic eects into account properly [2]. In order to
partially circumvent this problems, the authors introduced a threshold, roughly
corresponding to a molecule per reaction volume7, so that when the (continuous
valued) concentration falls below that level it is suddenly set to 0 [1].
In our model we stick to the original Kauman representation, but we relax
his non physical constraints that there be a maximum allowed length, and we
introduce the concentrations of the dierent types of molecules. The reactions
are described by stochastic dierential equations, which are solved by the well-
known Gillespie algorithm [8].
Moreover, we avoid three-molecular reactions because they are extremely rare.
So in our model we split the condensation reaction in two steps: rst the forma-
tion of a complex between the catalyst and a substrate and later the encounter
between the complex and the second substrate. The complex has a nite lifetime
so it may dissolve before condensation occurs.
The model is described in detail in section 3, where we describe also the problem
of the choice of the appropriate reaction graph.
Section 4 reports some results of our simulations, concerning in particular those
variables which aect the chance that an autocatalytic set is formed and main-
tained. In several cases one nds that there are key species with low numerosity
whose presence is crucial. This observation shows the usefulness of choosing a
7 The reaction volume was taken to be similar to that of a small bacterium, i.e. 1m3

4 Filisetti et al.
stochastic framework in order to follow the investigated processes.
In the nal section 5 indications for further work will be provided.
2 Critical Review of Previous Results
The main idea contained in the original work proposed in [12] is that the emer-
gence of autocatalytic sets of reactions is unavoidable when starting from a
mixture containing enough polypeptides. Considering polypeptides composed of
two monomers A and B and an initial population in which all polypeptides up
to length M are present, the total number of species is
SM =
MX
L=1
2L = 2M+1 2 : (1)
As anticipated in section 1, we consider two possible reactions: condensation
and cleavage. The total number of reactions building a specic polymer of length
L, 1 6 L 6M is then
RML;i =
MX
i=L+1
[2 2iL

] + (L 1) : (2)
It is possible to show [12] that the ratio between the total number of reactions
among polymers and the total number of species is equal to
RMtot
SM
=
MX
i=1
M i
2i
= [M 2] : (3)
Equation (3) shows that, although the total number of polymers increases
exponentially, the number of conceivable reactions increases even faster, leading
to the linear increase of their ratio. Therefore an autocatalytic set will certainly
form, no matter what the probability of catalysis is, provided that there are poly-
mers long enough, which in this case corresponds to having a sucient number
of dierent types of polymers in the system.
Equation (2) shows the number of reactions able to build a specic polymer of
length L, 1 6 L 6M . This number decreases as L approaches M , and there-
fore there are more ways to form short polymers than long polymers. However
one should also consider that for any given length L there are 2L dierent poly-
mers formed by a 2 letters alphabet; as consequence the number of reactions
which give rise to polymers of length L is
RML;tot = 2
L
"
MX
i=L+1
[2 2iL] + (L 1)
#
2L+1 : (4)

A stochastic model of autocatalytic reaction networks 5
the last term takes into account the case where the index of the sum is
i = 2L8). Therefore, although there are more ways to create a short specic
polymer than a longer one, the formation of long polymers is favoured.
The concentrations of the various species are taken into account in the models
[6, 1]. In both cases simulations were performed using a deterministic model, i.e.
a set of ordinary dierential equations, in which each equation describes the rate
of change of the concentration of a species.
Analyses were based on the connectivity ratio between edges and nodes, R,
of the reactions graph. If R < 0:5 the system is subcritical and no ACS is
formed, between 0:5 and 1 a critical zone is established and some structures
begin to form and R > 1 corresponds to the supracritical state, where cycles
appear. In particular they investigated the relations between the probability
that a randomly chosen molecule catalyses a randomly chosen reaction, and the
properties of the initial set of species, pointing out the importance of providing
a rich set of dierent species to allow the formation of cycles. Note that in this
model, like in that of Kauman, the initial food set is composed by all the species
up to length M, so it is impossible to distinguish the eects of the numerosity
from that of the maximum allowed length. If we call B the size of the alphabet
and L the maximum length of the initial ring disk, composed by all the possible
species, criticality is achieved if p  B2L. A a well stirred chemostat having
the same dimension of a bacterium is simulated by Bagley [1].
3 The Stochastic Model
3.1 The model
The aim of the model is not to provide a detailed description of a specic set of
reactions; rather, we want to focus our attention on the general characteristics
emerging from the interaction of a large number of interacting molecules. The
basic entities of this model are linear chains composed of a sequence of units
formed using a nite alphabet. The model is fairly general, and does not refer
to specic chemical classes; in particular, the basic units could represent sin-
gle elements, stable compounds or more complex domains. Given their essential
characteristic of allowing recursive assembly, let us call them (without losing
generality) \monomers", whereas the linear chains could be called \polymers"or
\species". Species will be used to denote either monomers or polymers.
Species can change in time, and they will be denoted at a given time by (s1; s2; : : : ; sn).
In the sequence of monomers the order matters, reading from left to right
(i.e. ABB is dierent from BBA) and the various species can have dierent
lengths. Any species can be present in multiple copies, and the number of ex-
emplars of the various species will be denoted by (x1; x2; : : : ; xn). In this work
the term \molecules" will be used to denote the number of exemplars of the
8 In this case one single reaction is enough to produce two dierent polymers of equal
length, so, the total amount of reactions have to be reduced by 2  22LL = 2L+1
reactions

6 Filisetti et al.
various species, either monomers or polymers. The reaction rates depend upon
concentrations but, working at a xed volume, we will sometimes refer to the
number of molecules.
As in the original model proposed by Kauman [12], we consider two possible
reactions, i.e. end-condensation and cleavage. Let us assume that both reactions
need a catalyst to occur and that the kinetic of the spontaneous reactions is
much slower than that of the catalysed reactions, hence spontaneous reactions
are neglected. We assume also that the rates of the reverse reactions are neg-
ligible with respect to that of the forward reactions. According to the number
and the length of the species present in the system the number of conceivable
reactions is
R =
NX
i=1
(L(si) 1) +N
2 : (5)
where L(si) is the length of the i-th species and N is the total number of
species present in the system (the sum refers to the number of cleavages and
the quadratic term to the number of condensations). Note that in any particular
simulation not all the possible reactions actually happen, and we assume that
there is a probability p that a randomly chosen species catalyses a randomly
chosen reaction. So in the present model there is no functional relationship be-
tween sequences (i.e. "chemical composition") of the catalysts and the reactions
they catalyse.
Concentrations change in time according to an asynchronous stochastic approach
based on the well-known Gillespie algorithm [8, 9] where at each step a reaction
is chosen (among all the possible ones) and physical time is computed. All reac-
tions occur in a well stirred chemostat and the concentration of each species is
assumed to be constant all over the space. The model behaviour thus depends
on the characteristics of the initial chemical population, on the composition of
the in
ow and on the value of the incoming
ux.
3.2 The Implementation
The model implementation is partitioned into two phases: initialisation and dy-
namics. The initialisation phase requires the creation of the initial population of
species. It is possible to build all the species having a length lower than or equal
to the initial maximum species length M (the so-called ring disks [6]), or to
generate a more sparse initial population, for example by removing species with
uniform probability till the desired initial number of species Nini (Nini < M)9
is reached, or by directly creating only Nini species uniformly choosing their
length l (l 2 [1;M ]) and composition. In the simulations described in this paper
all the species of length 1 are always created.
Once a species is created, its amount is initialised according to some distribu-
tion having the length of the species as a parameter; so it is possible to generate
9 In the following simulations we use initial ring disks and incoming
uxes without
the presence of ACS

A stochastic model of autocatalytic reaction networks 7
initial ring disks favouring short species or having a uniform distribution of
species lengths. The total number R of conceivable reactions is computed, and
each species will catalyse each reaction independently with probability p (there-
fore generating on average r = R
p catalysed reactions). The cutting point of
cleavages is chosen with uniform probability and it is possible to tune the rela-
tive fraction of cleavages and condensations.
In condensation reactions two substrates are selected and the product is created
binding them together. Note that, while cleavage is a bi-molecular reaction that
occurs instantaneously when the catalyst binds the substrate, condensation is
a three-molecular reaction and it is performed in two steps: in the rst step
the catalyst binds a substrate forming a temporary complex while in the second
step the complex binds the second substrate releasing product and catalyst. Note
that also the spontaneous degradation of the complex is considered. Hence the
system is composed of four possible kind of reactions (condensation being split
in three phases):
{ Cleavage: AB + C ! A+B + C
{ Condensation: (whole phenomenon: A+B + C ! AB + C)
 Complex formation: A+ C ! AC
 Complex dissociation: AC ! A+ C
 Final condensation: AC +B ! AB + C
The complex degradation and nal condensation do not require a specic
catalyst.
The outgoing
ux is simulated by means of a spontaneous decay of each specie
and complex, with a common kinetic constant. The incoming
ux is composed
by a xed number of species, each one having a denite concentration; molecules
enter one at a time at xed rate.
The internal dynamic must take into account the fact that new species can be
created; these new objects could react with the already present chemicals, and/or
catalyse new reactions. So, whenever a new species is created, N increases, and
the total number of conceivable reactions increases too; as consequence reac-
tions catalysed by the species already present in the systems have to be updated
according to the presence of the new species. If after a while the new specie
disappears, the new set of reactions is nevertheless kept in memory, to assure
consistency if that species reappears.
The model allows the emergence of competition and inhibition phenomena. The
former are related to the impossibility for a single catalyst to be involved in more
than one reaction at a time when it is forming an intermediate complex during
a condensation reaction, while the latter occur, for instance, when a component
of a reaction is consumed by other reactions thus decreasing the rate of the rst
reaction. In such a way the system can regulate its internal activity.
The adoption of an asynchronous updating framework creates a problem in iden-
tifying cycles. In order to analyse the structure of the system in fact it is useful
to represent species and reactions by means of a directed graph, where the nodes
are the chemicals, and node X is linked to node Y if the species corresponding to

8 Filisetti et al.
node X catalyses the formation of the specie corresponding to node Y . However,
because of the asynchronous updating of the algorithm, there can be several net-
work representations (a situation not encountered in the case of deterministic
approaches, where all the reactions occur at the same time, corresponding to
stable edges in the reaction graph representation).
Strictly speaking, on the graph representing the asynchronous framework at
each step there is only one edge (representing the single reaction occurred dur-
ing the last time interval). In order to avoid such a meaningless representation
two strategies are possible. The rst one involves the introduction of a temporal
window, so the graph contains all the reactions which have taken place in the
window. In order to neglect the in
uence of very rare reactions, i.e. those with
a very low rate but which might have occurred once in the past, each link has a
decay time: if the reaction does not occur within this time, the edge is removed
from the graph (\actual reaction graph" representation).
Another possible representation could be created by drawing a graph containing
all the possible reactions according to the species present in the system, with at
least one molecule, at a given time t. This representation could provide indica-
tions about the \nearest adjacent possible" [13].
In order to answer the question whether a cycle has been formed the rst rep-
resentation provides a more meaningful information; indeed in nature all the
reactions take place in parallel and placing them in a sequence is just an algo-
rithmic artefact.
4 Simulation Results
An ACS is a subgraph in which each node is directly or indirectly connected with
each other node belonging to the ACS, i.e. it is a strongly connected component
(SCC for short). This fact implies the presence of a path from a to b but also
from b to a, and as consequence a SCC contains at least one cycle.
In order to determine whether an autocatalytic set is present we use two dierent
methods related to the proprieties of the adjacency matrix A created by using
the \actual reaction graph" representation described above (the Aij elements is
equal to 1 when species j catalyses the formation of species i and 0 if it does
not).
According to the rst method one directly looks for the presence of strongly con-
nected components or of cycles in the reaction graph, while the second method is
based on the analysis of the eigenvalue of A with largest real part (ELRP), let us
denote it as 1. In fact, since A is real and non-negative, the Perron-Frebonius
theorem [14] assures that 1 is real and 1 > 0; it can also be shown that if
1 = 0 there are no cycles in the graph, and the presence of a cycle is associated
with 1 > 1 [10].
A rst group of simulations shows the dependence of the system behaviour on

10 Filisetti et al.
0 20 40 60 80 1000
200
400
600
800
Time
# of Mo
lecules
in SCC
s
0 20 40 60 80 1000
204
060
8010
0
Time
% of M
olecule
s in SC
Cs
Fig. 3. On the left panel the number of molecules belonging to the SCC is shown
whereas, on the right panel, the percentage of molecules (over the total number of
molecules) belonging to a SCC versus time is shown. All simulations are made using a
ring disk containing all the species up to length 4.
ACSs in the in
ux. So the existence of SCCs requires species which are gener-
ated in the reactor, and some of them are often present with a small number of
molecules. This fact has an important consequence: in fact, species characterised
by low number of molecules are subject to large
uctuations (see also gure 4)
which can in turn aect also species having greater numerosities, but that are
catalysed by rarer ones.
The results can be summarised saying that SCCs, and ACSs too, could be sub-
ject to heavy turbulence periods, and that simulations performed using a de-
terministic framework could not capture some important characteristics of the
phenomenon. Note that the in
uence of species having few molecules on ACS
maintenance could help in explaining some diculties associated in creating
ACSs in the laboratory.
0 20 40 60 80 1000
5
10
15
Time#
of Mol
s (NO IN
FLUX) in
SCCs
0 20 40 60 80 1000
204
060
8010
0
Time%
of Mol
s (NO IN
FLUX) in
SCCs
Fig. 4. On the left panel the number of molecules belonging to the SCC those are not
present in the in
ux is shown whereas, on the right panel, their percentage (over the
total number of molecules not present in the in
ux) belonging to a SCC versus time
is shown. All simulations are made using a ring disk containing all the species up to
length 4.

A stochastic model of autocatalytic reaction networks 11
5 Conclusions and Indications for Further Work
The model which we describe here allows a proper treatment of the stochas-
tic dynamics associated to collectively autocatalytic reaction sets, and therefore
provides a virtual laboratory for studying the eects of various features of the
reaction set, and of the values of the parameters. For example, it can be seen
that the chemical composition of the incoming
ux heavily in
uences the emer-
gence of ACSs, while the composition initially present in the reactor has a minor
eect.
The use of a stochastic approach is well grounded a priori, but it is also conrmed
by the results of the simulations. Indeed one of the most interesting outcomes
is the fact that, in most cases, the existence of the ACS critically depends upon
that of some rare molecules, which may disappear due to
uctuations. The model
might be used in the future to analyse possible strategies to get rid of these in-
stabilities.
Many improvements are of course possible, like for example that of relating the
structure (i.e. the chemical composition) of a catalyst to the kind of reactions
it can in
uence. Moreover, the model presented here has essentially an informa-
tional nature, and energy is not considered explicitly. Since energy plays a major
role in nature, a priority in future developments is that of integrating it in our
general picture of evolving reaction sets.
Acknowledgments. This work has been partially supported by the Fondazione
di Venezia, http://www.fondazionedivenezia.it, (DICE project). Interesting dis-
cussions with Davide De Lucrezia are gratefully acknowledged.
References
1. R.J.Bagley and J.D.Farmer. Spontaneous emergence of a metabolism. Articial
Life II: 93-140. Santa Fe Institute Studies in the Sciences of Complexity, X, 1991.
2. J. M. Bower and H. Bolouri. Computational Modeling of Genetic and Biochemical
Networks (Computational Molecular Biology). The MIT Press, 2004.
3. T. Carletti, R. Serra, I. Poli, M. Villani, and A. Filisetti. Sucient conditions for
emergent synchronization in protocell models. J Theor Biol, 254(4):741{751, 2008
Oct 21.
4. F. J. Dyson. Origins of life. Cambridge, England, 1985.
5. M. Eigen and P. Shuster. The hypercycle: A principle of natural self-organization.
Springer, Berlin, 1979.
6. J. D. Farmer, S. A. Kauman, and N. H. Packard. Autocatalytic replication of
polymers. Physica D, 22(2):50{67, 1986.
7. A. Filisetti, R. Serra, T. Carletti, M. Villani, and I. Poli. Synchronization phenomena
in protocell models. Biophysical Reviews and Letters (BRL), 3(1/2):325{342, April
2008.
8. D. T. Gillespie. Exact stochastic simulation of coupled chemical reactions. The
Journal of Physical Chemistry, 81(25):2340{2361, 1977.

12 Filisetti et al.
9. D. T. Gillespie. Stochastic simulation of chemical kinetics. Annual Review of Phys-
ical Chemistry, 58(1):35{55, 2007.
10. S. Jain and S. Krishna. A model for the emergence of cooperation, interdependence,
and structure in evolving networks. Proc Natl Acad Sci U S A, 98(2):543{547, 2001
Jan 16.
11. K. Kaneko. Life: An Introduction to Complex Systems Biology (Understanding
Complex Systems). Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2006.
12. S. A. Kauman. Autocatalytic sets of proteins. J Theor Biol, 119(1):1{24, 1986.
13. S. A. Kauman. Reinventing the sacred : a new view of science, reason and religion
/ Stuart A. Kauman. Basic Books, New York :, 2008.
14. H. Lutkepohl. Handbook of Matrices. New York: John wiley & Sons, xvi + 304
pp., ISBN 0-471-96688-6, 1996.
15. A. Munteanu, C. Attolini, S. Rasmussen, H. Ziock, and R. Sole. Generic darwinian
selection in protocell assemblies. DOI: SFI-WP 06-09-032, SFI Working Papers,
Santa Fe Institute, 2006.
16. S. Rasmussen, L. Chen, D. Deamer, D. C.Krakauer, N. H. Packard, P. F. Stadler,
and M. A. Bedau. Transitions from nonliving to living matter. Science, 303, 963-965,
2004.
17. R. Serra, T. Carletti, and I. Poli. Syncronization phenomena in surface-reaction
models of protocells. Articial Life 13: 1-16, 2007.
18. J. W. Szostak, D. P. Bartel, and P. L. Luisi. Synthesizing life. Nature,
409(6818):387{90, Jan 2001.
19. N. Takeuchi and P. Hogeweg. Multilevel selection in models of prebiotic evolution
ii: A direct comparison of compartmentalization and spatial self-organization. PLoS
Comput Biol, 5(10):e1000542, 10 2009.