The influence of the residence time on the
dynamics of catalytic reaction networks
Alessandro FILISETTI a,1, Roberto SERRA b,a, Marco VILLANI b,a,
Alex GRAUDENZI a, Rudolf M. FÜCHSLIN a,c, and Irene POLI a,d
a European Centre for Living Technology, S.Marco 2940, 30124 - Venice (IT),
http://www.ecltech.org
b Dipartimento di Scienze Sociali, Cognitive e Quantitative
Università di Modena e Reggio Emilia, via Allegri 9, 42100 Reggio Emilia, Italy
c Artificial Intelligence Lab Univ. Zürich
Andreasstr. 15, CH-8050 Zürich, Switzerland
Office: +41 (0)44 635 45 9, Mobile: +41 (0)79 232 74 36
d Dipartimento di Statistica, Università Ca’ Foscari,
San Giobbe - Cannaregio 873, 30121 Venezia, Italy
Abstract. Although autocatalytic networks are common in nature, it is very diffi-
cult to reproduce them in laboratory. Since there are several models in literature
describing a phase transition to an autocatalytic set once that a certain degree of
heterogeneity in the composition of the system is reached, it is interesting to un-
derstand why it is so difficult to observe such a phenomenon in the laboratory. For
this reason, we here present a model designed for the study of that systems taking
into account the stochastic nature of the dynamics of interacting molecules. In par-
ticular, the analysis is focused on the emergence of autocatalytic sets in accordance
with different residence times and influx compositions.
Keywords. Complex systems biology, Catalytic reaction networks, Autocatalytic
sets of molecules
Introduction
The discovery and the description of the generic properties of those systems characterised
by the emergence of sets of autocatalytic networks turn out to be of central importance
in the investigation on the origin of life, as well as in the theoretical design of artificial
protocells [3,8,18,19,22].
To this end, in the last decades many different models, with distinct features and objec-
tives, have been developed, among which the well-known works by Dyson [5], Eigen
and Schuster [6], Kauffman [15], Jain and Khrishna [13] and Kaneko [14].
The models differ with regard to the hypotheses concerning the physical nature of the au-
tocatalytic networks in the real world, describing either gene-first or protein-first scenar-
ios [4,24]. On the other hand, a wide consensus has recently emerged about the necessity
of a container that separates the set of collectively set of molecules from the environment
1Corresponding Author, E-mail: alessandro.filisetti@ecltech.org.
Neural Nets WIRN10
B. Apolloni et al. (Eds.)
IOS Press, 2011
© 2011 The authors and IOS Press. All rights reserved.
doi:10.3233/978-1-60750-692-8-243
243

and whose growth and division rates are coupled and synchronised with the dynamics of
the collectively self-replicating sets of molecules [3,10,20,21,14,18,22].
We here introduce a model that is an evolution of the model originally proposed by
Kauffman [15]. The original model describes a system that is composed by linear chains
of monomers symbolised with letters taken from a finite alphabet. Only two possible
reactions are allowed, namely condensation, in which two monomers or polymers are
joined together, and cleavage, in which a polymer is divided in two parts. All the reac-
tions occur at a negligible rate unless they are catalysed by some molecule belonging
to the system. On the basis of combinatorial arguments regarding the graph of all the
conceivable reactions, Kauffman stated that, whatever is the probability for a reaction
to be catalysed by a molecule present in the system, increasing the maximum length of
the initial set of molecules, the number of conceivable reactions grows faster than the
number of molecules, hence, at least an autocatalytic set (ACS from now on) is going to
emerge.
Nevertheless, since the original Kauffman model neither take into account the actual
concentration of the molecules, nor even possible selection phenomena involved in the
process, Farmer and others [7,1] designed a model in which the system is described by a
set of differential equations, thus introducing a form of dynamics, and modelled within a
chemostat, which may provide a sort of selective pressure. The results of the simulations
showed the existence of a range of structural parameters according to which at least an
ACS actually emerges, hence supporting the hypothesis of Kauffman. In particular they
recognised that the emergence of an ACS is correlated with the average connectivity of
the reactions graph: below a value of 0.5 nothing happens, whereas for values larger than
1 different ACSs emerge and in the middle they observed a critical region in which dif-
ferent clusters of connected molecules begin to form. Thus, considering that the average
connectivity is a function of the probability for a reaction to be catalysed, for any given
probability there is always a critical size of the initial set of molecules so that strongly
connected components begin to form.
In spite of the deterministic description proposed by Farmer et al [7,1], we introduce
an asynchronous stochastic model, implemented by means of the well-known Gillespie
algorithm [11], in order to consider the actual amount of each molecule and the influ-
ence that such molecules present in a few copies may have on the dynamics. Other chief
novelties are the removal of the constraints regarding the maximum allowed length of
the molecules and the distinction between bi-molecular (cleavage) and three-molecular
(condensation) reactions. The condensation process is split in two steps, the formation
of a complex, characterised by a finite lifetime, involving the catalyst and a substrate and
the subsequent encounter between the complex and a second substrate.
In section 1 the model is described, alongside with the considerations about the choice of
the appropriate reaction graph representation. In section 2 the results obtained simulating
a simple substrate-enzyme kinetic, both with a deterministic and a stochastic description,
are shown. In section 3 the results of the computational simulations in which the resi-
dence time is changed are presented, while in the final section the conclusions and the
indications for further works are provided.
A. Filisetti et al. / The Influence of the Residence Time244

1. Description of the model
For an exhaustive description of the model, the reader is referred to [9]. Here we will
only outline its main features in order to provide a clear picture of the system.
The foremost aim of the model is the detection and the characterisation of the general
properties characterising the dynamics emerging from a large number of molecules that
interact through catalysed reactions. The basic entities of the system are chains of ele-
ments represented with the letters of a finite alphabet. Since the level of abstraction is
high and the model general, there is not specific correspondence among the letters and
any particular real chemical class: on the contrary, basic units may represent single el-
ements, stable compounds or more complex domains. Since one of the key properties
is the possibility of recursive assembly, we will define them as “monomers”, while the
ordered linear chains will be defined as “polymers” (e.g. AAB is different from BAA).
“Species” will denote both monomers and polymers.
The set of species is not fixed in time and is denoted at a given time by the list
S = (s1, s2, ..., sn). Any species can be present in the system at a certain time in multi-
ple copies, defined as “molecules” and whose number is given by X = (x1, x2, ..., xn).
We consider only two kinds of reactions, as in the original Kauffman model, condensa-
tion and cleavage. We assume that both the condensation and the cleavage reactions need
a catalyst (i.e. a species) to occur, while the spontaneous reactions are neglected. Simi-
larly, we assume that there are forward reactions (of both the kinds), while the backward
reactions are in principle ignored.
Each species is chosen to be a catalyst of a particular reaction according to an indepen-
dent probability p; therefore neither all the conceivable reactions will actually occur, nor
will all the species will be catalysts. Notice also that a species can be a catalyst for more
than one reaction and that a reaction can be catalysed by more than one species.
According to the number and the length of the species present in the system, the number
of conceivable reactions is:
R =
N
∑
i=1
(L(si)− 1) + N
2 . (1)
where L(si) is the length of the i − th species and N is the total number of species
present in the system (the first term refers to the number of cleavages and the second to
the number of condensations).
On the basis of R, each species will catalyse r = R × p reactions (on the average). No
functional relationship between sequences (i.e. “ chemical composition”) of the catalysts,
substrates and the reactions is considered in the model.
Numerosities are introduced within the system according to an asynchronous stochastic
approach based on the Gillespie algorithm [11,12]: at each step a reaction is chosen at
random among all the possible ones and its physical time is computed. All reactions
occur in a well stirred chemostat and the concentration of each species is assumed to be
constant in space.
1.1. The Implementation
The implementation of the model entails two distinct phases: initialisation and dynamics.
In the initialisation phase the initial population of species is created (the so-called firing
A. Filisetti et al. / The Influence of the Residence Time 245

disk [7]). It is possible to create a) all the species whose length is equal or lower than
the initial maximum species length M , b) a more sparse initial population, for instance
by the removal of a certain number of species from the set of all the possible species, c)
directly only a number Nini species choosing their length and composition. In the simu-
lations that we will discuss in section 3 we decided to create in any cases all the possible
monomers.
Once that all the species of the firing disk have been created, their amount is initialised in
accordance with some frequency distribution related to their length. It is possible to gen-
erate firing disks in which a) short species are favoured, b) there is a uniform distribution
of species lengths. Each species will be chosen to be a catalyst of a random reaction with
uniform probability p and it is possible to favour either cleavages or condensations.
Given two species, for instance A and B, standing for substrates and a third species C,
standing for a catalyst, the scheme of the two possible reactions is the following:
• Cleavage: AB + C → A + B + C
• Condensation: (whole phenomenon: A + B + C → AB + C)
∗ Complex formation: A + C → C.A
∗ Complex dissociation: C.A → A + C
∗ Final condensation: C.A + B → AB + C
Condensation is split in two distinct phases: a) the complex formation/dissociation,
in which the catalyst binds a substrate (A) forming a temporary complex (C.A) (i.e. for-
ward reaction) or vice-versa (i.e. backward reaction), b) the final condensation in which
the complex binds the second substrate (B), releasing the product (AB) and the catalyst.
Notice that the cutting point of cleavages is chosen with uniform probability.
To simulate the outgoing flux, a spontaneous decay is associated to each species, with a
common kinetic constant. Conversely, the ingoing flux is composed by species belonging
to a sub-set of the initial firing disk, each one with a specific concentration.
New species can be created and become a part of the internal dynamics, reacting with
the species already present in the system and/or catalysing novel reactions: when a new
species is created, N increases, together with the number of possible reactions, ensuring
the coherence of the already existing reactions scheme.
To analyse the structure of the system, an effective instrument is the directed graph in
which the nodes are the species and there is a link from X to Y if X is a catalyst of
a reaction (cleavage or condensation) in which Y is a product. Nevertheless, since the
updating mechanism is asynchronous, it is difficult to unambiguously identify cycles. In
order to overcome this setback, we introduce a specific decay time associated to each
reaction in order to neglect the influence of very rare reactions: if the reaction does not
occur within this time, the edge is removed from the graph (this is defined as the “actual
reaction graph” representation).
Notice also that the system tends to regulate its internal dynamics through emergent phe-
nomena of competition and inhibition, respectively due to the exclusivity of the catalysts
involved in specific reactions and on the possible consumption of some species involved
in more than one reaction.
A. Filisetti et al. / The Influence of the Residence Time246

2. Simple substrates-enzyme kinetic simulation
Although the stochastic simulator has been designed to simulate a large set of interacting
molecules it is effective in describing the dynamics of simple systems as well. For this
reason a simple substrates-enzyme kinetic system has been simulated, using a determin-
istic description through a set of differential equations, numerically integrated, to com-
pare with the results obtained with the stochastic simulator, using the same parameters.
Let molecule E stands for a general enzyme, which catalyses a reaction, in a fixed volume
with thermal equilibrium, wherein the substrates S1 and S2 are converted into product
P . We consider also the reverse reaction in which P is converted in S1 and S2, always
by means of the catalyst E. Therefore, the complete reaction scheme is the following:
• S1 + S2 + E  P + E
An intermediate step of the reaction is the creation of a complex C formed by means
of the collision between one of the two substrates (i.e. S1) with the enzyme E. Hence,
the above reaction scheme can be split in:
a) S1 + E  C + S2 → P + E
b) P + E → E + S1 + S2
The temporary complex C can be involved in two different reactions: a) binding
with the other substrate S2 and releasing the final product P and the enzyme E, b)
disassociating into enzyme E and substrate S1. Therefore, four distinct reactions are
present in the system:
1. S1 + E → C with rate k1,
2. C → S1 + E with rate k−1,
3. C + S2 → P + E1 with rate k2,
4. E + P → E + S1 + S2 with rate k3.
The system is then described by the following set of ordinary differential equations:
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎩
d[S
1
]
dt = k−1[C]− k1[S1][E] + k4[E][P ]
d[S
2
]
dt = −k2[C][S2] + k3[E2][P ]
d[E]
dt = k−1[C] + k2[C][S2]− k1[S1][E]
d[C]
dt = k1[S1][E]− k−1[C]− k2[C][S2]
d[P ]
dt = k2[C][S2]− k3[E][P ]
(2)
Notice that if we set constant k3 = 0 the reverse reaction is neglected: in regard
to this case, the graphs in fig. 1 show that the behaviour of the stochastic simulation,
fig. 1(b), well approximates that provided by the numerical integration of eqs. 2, fig. 1(a).
Results regarding the presence of the reverse reaction are, instead, shown in fig. 2.
As one can see, the results of the stochastic simulation compared with those obtained by
means of the deterministic approach are indeed similar.
A. Filisetti et al. / The Influence of the Residence Time 247

0 10 20 30 40 500
1
2
3
4
5x 10
−7
Time
Co
nc
en
tra
tio
n


[S1]
[S2]
[P]
[E]
(a) 0 10 20 30 40 50
0
1
2
3
4
5x 10
−7
Times
Co
nc
en
tra
tio
n

[S1]
[S2]
[E]
[P]
(b)
Figure 1. On the left panel the deterministic description of a simple substrates-enzyme kinetic system with-
out reverse reaction is shown while on the right panel the behaviour obtained by the stochastic simulator is
shown. On the X axis time is represented whereas, on the Y axis, the concentration is represented (Parameters:
S
1
(0) = 5e − 7, S
2
(0) = 5e − 7, E(0) = 2e − 7, C(0) = 0, P (0) = 0, k
1
= 1e6, k
−1
= 1e − 4,
k
2
= 6e5, K
3
= 0).
0 10 20 30 40 500
1
2
3
4
5x 10
−7
Time
Co
nc
en
tra
tio
n


[S1]
[S2]
[P]
[E]
(a) 0 10 20 30 40 50
0
1
2
3
4
5x 10
−7
Time
Co
nc
en
tra
tio
n


[S1]
[S2]
[E]
[P]
(b)
Figure 2. On the left panel the deterministic description of a simple substrates-enzyme kinetic system with
reverse reaction is shown while on the right panel the behaviour obtained using the stochastic simulator is
shown. On the X axis time is represented whereas, on the Y axis, the concentration is represented (Parameters:
S
1
(0) = 5e − 7, S
2
(0) = 5e − 7, E(0) = 2e − 7, C(0) = 0, P (0) = 0, K
1
= 1e6, K
−1
= 1e − 4,
K
2
= 6e5, K
3
= 1e6).
3. Results of the simulations of systems with different residence times
In our previous work [9] the investigations have been focused on the behaviour of the
system according to different influx compositions in terms of species diversity. The main
result is that in correspondence of an incoming flux characterised by a raising diversity,
the probability that the reaction graph contains autocatalytic sets increases. It is impor-
tant to remark that we are referring to the structural proprieties of the reaction graph and
that such proprieties do not take into account the possible consequences of the presence
of an ACS on the system dynamics. Following [6,1], a set of autocatalytic molecules
should produce a significative departures from the equilibrium concentration in order to
yield robustness and evolvability.
In this paper we address our analysis on the contribute provided by the residence time
of the molecules within the reactor. In our previous simulations [9] we considered an
average residence time of 10 seconds for the molecules in the reactor with an incoming
flux of 100 mols/sec, whereas in this work the results obtained in [9] have been compared
with those obtained in accordance with a residence time increased to 100 seconds and an
influx of 1000 mols/sec.
The idea is that increasing the residence time and the concentration of the molecules
composing the incoming flux, the chance for the molecules to collide with each other
increases as well. It is interesting to investigate whether the enhanced activity of the
A. Filisetti et al. / The Influence of the Residence Time248

system provides the same results or some possibly new emergent behaviour.
The graphs in fig. 3 show the cumulative number of species, i.e. the new species, in
function of time. The simulations characterised by a residence time of 10 seconds clearly
show, fig. 3(a), 3(b) and 3(c), an increment in new species creation according to a rais-
ing diversity, whereas the systems with a longer residence time, fig. 3(d), 3(e) and 3(f),
corresponding to an influx compose of all the molecules up to length 2, 3 and 4, does not
show particular differences in new species production.
     





































Figure 3. The graphs represent the average cumulative number of species, i.e. the new species, in function
of time, over ten different simulations (the error bars represent the standard deviation). On the top panel,
respectively figures (a), (b) and (c), the results of the simulations regarding the system characterised by an
average residence time of 10 seconds and an incoming flux composed of 100 mols/sec are represented. On
the bottom panel, respectively figures (d), (e) and (f), the results of the simulations concerning the system
characterised by an average residence time of 100 seconds and an incoming flux composed of 1000 mols/sec
are represented.
The same behaviour is observed also looking at both the total number and the max-
imum lengths of the living species (not shown here). It seems that increasing the resi-
dence time the influence of the diversity in the incoming flux does not affect the over-
all behaviour of the system. We may suppose that when the species have a shorter time
to react with each other the only chance to have a self maintaining system is feeding it
with a rich set of molecules. On the contrary, a longer residence time seems to confer an
higher degree of stability to the system, reason that the system survives also in the case
in which only species up to length 2 are provided.
Another interesting behaviour is shown in the graphs of fig. 4. Although the system with
a short residence time always shows different behaviours in accordance with the influx
diversity, and considering a longer residence time the behaviours of the different influxes
are essentially the same, the performance associated feeding it with all the molecules up
to length 4 seems to be higher if compare to that of systems with a longer residence time.
It may be that increasing the residence time, and the species concentration, the system
is, on the one hand, able to self-maintaining feeding it with less nutrients but, on the
other hand, the enhanced activity of the system provides not only useful reactions but
also those that tends to destroy the autocatalytic sets.
A. Filisetti et al. / The Influence of the Residence Time 249

Acknowledgements
This work has been partially supported by the Fondazione di Venezia, http://ww .
(DICE project). Interesting discussions with Davide De Lucrezia
edged.
References
[1] R.J.Bagley and J.D.Farmer. Spontaneous emergence of a metabolism. Artificial 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 (Compu-
tational Molecular Biology). The MIT Press, 2004.
[3] T. Carletti, R. Serra, I. Poli, M. Villani, and A. Filisetti. Sufficient conditions for emergent synchroniza-
tion in protocell models. J Theor Biol, 254(4):741–751, 2008 Oct 21.
[4] D. De Lucrezia, F. Anella, and C. Chiarabelli. Question 5: on the chemical reality of the rna world. Orig
Life Evol Biosph, 37(4-5):379–385, 2007.
[5] F. J. Dyson. Origins of life. Cambridge, England, 1985.
[6] M. Eigen and P. Shuster. The hypercycle: A principle of natural self-organization. Springer, Berlin,
1979.
[7] J. D. Farmer, S. A. Kauffman, and N. H. Packard. Autocatalytic replication of polymers. Physica D,
22(2):50–67, 1986.
[8] 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.
[9] A. Filisetti, R. Serra, M. Villani, R. M. Füchslin, N. H. Packard, S. A. Kauffman, and I. Poli. A stochastic
model of autocatalytic reaction networks. Submitted to European Conference on Complex Systems,
2010.
[10] A. Filisetti, R. Serra, T. Carletti, M. Villani, and I. Poli. Non linear reaction models of protocell, syn-
chronization and chaos. European Physical Journal B (in press), 2010.
[11] D. T. Gillespie. Exact stochastic simulation of coupled chemical reactions. The Journal of Physical
Chemistry, 81(25):2340–2361, 1977.
[12] D. T. Gillespie. Stochastic simulation of chemical kinetics. Annual Review of Physical Chemistry,
58(1):35–55, 2007.
[13] 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.
[14] K. Kaneko. Life: An Introduction to Complex Systems Biology (Understanding Complex Systems).
Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2006.
[15] S. A. Kauffman. Autocatalytic sets of proteins. J Theor Biol, 119(1):1–24, 1986.
[16] S. A. Kauffman. Reinventing the sacred : a new view of science, reason and religion / Stuart A. Kauff-
man. Basic Books, New York :, 2008.
[17] H. Lutkepohl. Handbook of Matrices. New York: John wiley & Sons, xvi + 304 pp., ISBN 0-471-96688-
6, 1996.
[18] 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.
[19] R. Serra, T. Carletti, and I. Poli. Syncronization phenomena in surface-reaction models of protocells.
Artificial Life 13: 1-16, 2007.
[20] R. V. Solé, J. Macía, H. Fellermann, A. Munteanu, J. Sardanyés, and S. Valverde. Models of protocell
replication. Protocells: bridging living and non-living matter, pages 213–231, 2008.
[21] R. V. Solé, A. Munteanu, C. Rodriguez-Caso, and J. Macía. Synthetic protocell biology: from reproduc-
tion to computation. Philos Trans R Soc Lond B Biol Sci, 362(1486):1727–39, Oct 2007.
[22] J. W. Szostak, D. P. Bartel, and P. L. Luisi. Synthesizing life. Nature, 409(6818):387–90, Jan 2001.
[23] N. Takeuchi and P. Hogeweg. Multilevel selection in models of prebiotic evolution ii: A direct com-
parison of compartmentalization and spatial self-organization. PLoS Comput Biol, 5(10):e1000542, 10
2009.
[24] G. Wächtershäuser. From volcanic origins of chemoautotrophic life to bacteria, archaea and eukarya.
Philos Trans R Soc Lond B Biol Sci, 361(1474):1787–806; discussion 1806–8, Oct 2006.
fondazionedivenezia.it,
are gratefully acknowl
w.
A. Filisetti et al. / The Influence of the Residence Time 251