Dynamical Stability of Autocatalytic Sets
Rudolf M. Fu¨chslin1,3, Alessandro Filisetti1, Roberto Serra2,1, Marco Villani2,1,
Davide DeLucrezia1, Norman Packard4,1, Stuart A. Kauffman5,1 and Irene Poli6,1
1 European Centre for Living Technology
Calle del Clero 2940, 30124 Venice, Italy
2 Dipartimento di Scienze Sociali, Cognitive e Quantitative Universita di Modena e Reggio Emilia,
via Allegri 9, 42100 Reggio Emilia, Italy
3 Artificial Intelligence Lab Univ. Zu¨rich
Andreasstr. 15, CH-8050 Zu¨rich, Switzerland
4 ProtoLife Inc
57 Post St. Suite 513, San Francisco, CA 94104
5 Departments of Biochemistry and Mathematics University of Vermont,
Burlington, VT 05405
6 Dipartimento di Statistica, Universita Ca’ Foscari, San Giobbe - Cannaregio 873,
30121 Venezia, Italy fuchslin@ifi.uzh.ch
Abstract
Theoretical investigations of autocatalytic sets rendered the
occurrence of self-sustaining sets of molecules to be a generic
property of random reaction networks. This stands in some
contrast with the experimental difficulty to actually find such
systems. Nature proves that autocatalytic sets can exist (and
in a robust manner), in fact they are ubiquitous, though usu-
ally operating under finely tuned conditions. In this work, we
argue that the usual approach, which is based on the study of
static properties of reaction graphs has to be complemented
with a dynamic perspective in order to avoid overestimation
of the probability of getting autocatalytic sets. Especially un-
der the, from the experimental point of view, important flow
reactor conditions, it is not sufficient just to have a pathway
generating a given type of molecules. The respective pro-
cess has also to happen with a sufficient rate in order to com-
pensate the outflow. Reaction rates are therefore of crucial
importance. Furthermore, processes such as cleavage are on
one hand advantageous for the system, because they enhance
the molecular variability, and therefore the potential for catal-
ysis. On the other hand, cleavage may also act in an inhibit-
ing manner by the destruction of vital components: therefore,
an optimal balance between ligation and cleavage has to be
found. If energy is included as a limiting resource, the con-
centration profiles of the components of autocatalytic sets are
altered in a manner that renders a certain range for the energy
supply rate as optimal for the realization of robust autocat-
alytic sets. Such a supply rate has to be chosen in order to get
a detectable signal in an experiment.
The results presented are based on a theoretical model and
obtained by numerical integration of systems of ODE. This
limits the number of involved molecular species which im-
plies that the quantitative findings of this work may have no
direct relevance for experimental situations. We claim, how-
ever, that the qualitative findings generalize (even in a more
pronounced manner) to systems of truly combinatorial size.
Keywords: Autocatalytic sets, emergence, origin of life
Introduction
In recent years, autocatalytic sets (ACS) have attracted inter-
est from many different research directions. Probably most
prominent are thereby investigations concerning the origin
of life, but ACS proved to be a concept also of value e.g.
for the study of transitions in general (non-chemical) sys-
tems of interacting production processes (6), including the
generation of knowledge.
Informally, the fundamental question with respect to
chemical reaction networks is whether or not a given set
of different, potentially catalytic molecules immersed into
a suitable environment (most often some type of flow reac-
tor) and provided with a sufficient supply of food or building
blocks is able of maintaining the concentration of its mem-
bers via mutual catalysis. The conditions under which such a
self-maintaining or autocatalytic set can be expected to ap-
pear with sufficiently high probability are then those to be
mimicked in an experiment e.g. concerned with the emer-
gence of protolife.
Based on different models of catalytic networks, there
is broad literature (10; 12; 7) on the detection of ACS ex-
ploiting combinatorial properties of reaction graphs. In (7)
a polynomial-time algorithm for the detection of a specific
type of ACS has been presented. Hordijk and Steel applied
this algorithm to a model by Kauffman (9). By analyzing
large numbers of randomly chosen networks, they corrobo-
rated a conclusion (to be discussed in detail below) which
Kauffman derived from combinatorial reasoning, namely
that in sufficiently diverse populations of potentially cat-
alytic chain molecules ACS will emerge. Thereby, ACS
will form independent of how sparse catalytic activity is dis-
tributed in the combinatorial variety of molecules, as long
as this variety is big enough (usually limited by a maximal

sequence length). Stated differently, given a certain variety
of potentially catalytic molecules, there is always a thresh-
old for the probability of catalytic activity such that above
that threshold, ACS can be expected to emerge with high
probability.
Despite some criticism (11) (for a discussion of Lifson’s
arguments, see (15)) and the fact that more detailled mod-
els of catalysis may modify some results presented in (9),
the main conclusions seem to generalize in one or the other
form to a broad variety of models. The obvious question to
ask then is, why ACS are not regularly discovered in the lab-
oratory. In (3), three possible answers were discussed. The
first one (sometimes preferred by experimentalists) claims
that the simplifications used in the formulation of the mod-
els on one hand make them tractable by analytical and/or
computational means but on the other hand renders them
unrealistic. The second answer (favored by some theorists)
says that the basic statements derived from simplified mod-
els are also valid if the details of the physical and chemi-
cal world were considered, but that the threshold necessary
for the emergence of ACS never has been reached. Finally,
the third position (and also the one advocated in (3) and in
this work) highlights the fact that in investigations purely
based on the properties of reaction graphs, dynamical and
stochastic aspects are not considered. For some models, this
is not necessary because their dynamics is basically (at least
piecewise) determined by linear operators (8). But for most
models based on general reaction graphs, graph-theoretical
methods may identify ACS which are only transient in the
sense that the chemical dynamics eventually leads to a col-
lapse of the ACS. This holds especially under flow reactor
conditions, where e.g. a catalyst needs not only to be pro-
duced via some reaction path, but also at a sufficient rate in
order to compensate for loss by outflow. Graph-theoretical
means are able to identify whether or not a reaction path is
present in a given network, but not wether the dynamics es-
tablishes a non-trivial stationary ACS (In fact, one should
speak of ACS exhibiting stationary or limit cycle behavior,
but in practice one observes most models to yield almost
exclusively stationary solutions. For a discussion, see e.g.
(14)). In an experiment, however, it may be difficult to ob-
serve transient ACS, first because they may only be active
during a very short period of time and second because their
emergence may be highly susceptible to initial conditions.
Stationary ACS which are able to produce a permanent de-
viation of some molecular concentrations from those one ex-
pects to result from the inflow and some non-catalytic back-
ground reactions offer a higher potential for being observ-
able in a reproducible manner (13).
Whereas in (3) the emphasis has been put on the investiga-
tion of the influence of stochastic fluctuations on the emer-
gence and dynamics of ACS, this paper is concerned with
the study of the influence of various parameters on the ob-
servabilty of stationary ACS. In section we discuss two dif-
ferent approaches for the definition of an ACS and motivate
the choice being taken for the investigations in this work. In
section , we briefly review the original model by Kauffman
(9) and present our implementation as a system of coupled
ODEs. In section , we show that the emergence of station-
ary ACS depends critically on the choice of parameters. Not
only the density of catalysts in the set of available molecules
is of importance, but also the relative occurrence of different
processes, say they ratio between ligation and cleavage pro-
cesses. We further study a derivative of the original model
that takes energy considerations into account, means the dif-
ferent reactions compete for a, with a constant rate renewed,
energy resource. We close with a discussion of the relevance
of our results for experimental setups.
Autocatalytic Sets
We compare two different approaches for the definition of
an autocatalytic set. The first one is especially appropriate
for the study of reaction graphs and thoroughly discussed
and formalized e.g. in (7). The second one, e.g. dis-
cussed in (13) takes into account the dynamics of the system
but is less formal. Bagley and Farmer define an “autocat-
alytic metabolism” as a coupled set of reactions which lead
to permanent concentrationswhich are significantly depart-
ing from the values one would obtain without catalysis. As
they point out themselves, this definition is to some extent
problematic, because what one regards as significant may
depend on the experimental means. However, we will use
a similar approach, because only those systems delivering a
measurable deviation (both with respect to quantities as well
as time) from some equilibrium distribution are of experi-
mental interest. In order to highlight the difference between
the two approaches, we briefly review the graph theoretical
definition used by Hordijk and Steel and show that an ACS
identified with their method needs not necessarily to be ob-
servable.
In (7) the main focus is laid on so called “reflexively au-
tocatalytic and F -generated reaction systems (RAF)” which,
informally, are those sets of reactions for which it holds that
a) each reaction is catalyzed by a molecule being part of the
reaction system and b) all reactants can be generated from
a food set F by iterative applications of the reactions in the
RAF.
In order to formalize the notion of a RAF, some defini-
tions are required (We adhere closely to (7)):
1. M denotes the set of molecules, e.g. a combinatorial va-
riety of chain molecules.
2. A reaction r is an ordered pair (A,B) with A,B ⊂ M
and A ∩ B = ∅. Thereby, A represents the reactants
and B the products of a reaction n1a1 + . . . + nkak →
n′1b1 + . . .+ n′lbl, ai ∈ A, bj ∈ B and the n, n′ represent
positive integers.

setting.
The Basic Model
In (9) the properties of sets of potentially catalytic di-block
copolymers were investigated. Thereby, it was assumed
1. Polymers consist of two different types of monomers A
and B.
2. There are two types of catalyzed reactions, namely liga-
tion and cleavage.
3. The probability for a polymer pc to catalyze a ligation
p1 + p2
pc−→ p1p2 or a cleavage p1p2
pc−→ p1 + p2 is
given by a probability r. This means that in the course
of the initialization of the network, catalytic processes are
set once and kept during an individual run.
This setting, basically a random reaction system, doesn’t
make any specific “helpful” assumptions supporting the
emergence or existence of an ACS, and nevertheless, strong
evidence was given that such a system should eventually
contain ACS, given only a sufficiently large variety of dif-
ferent polymers being included in the system (In case of
block polymers, this can be achieved simply by allowing se-
quences of length up to a critical Lc).
Several implementations of random graph models using
ODEs have been studied, see e.g. (1; 13). In this work, the
dynamics of the system is given by:
dpi
dt
= ki,in − dpi (2)
+ kj,k,L
∑
j,k,m
L(pj , pk, pi, pm)pjpkpm
− ki,j,L
∑
j,k,m
L(pi, pj , pk, pm)pipjpm
− kj,i,L
∑
j,k,m
L(pj , pi, pk, pm)pjpipm
+ kC
∑
j,k,m
C(pi, pj , pk, pm)pkpm
+ kC
∑
j,k,m
C(pj , pi, pk, pm)pkpm
− kC
∑
j,k,m
C(pj , pk, pi, pm)pipm.
Thereby, pi represents a polymer with given sequence com-
posed of two types of monomers A,B. The rate of in-
flux ki,in under flow reactor conditions is set to one for the
monomers pi = A,B and zero for all other sequences. Out-
flow is determined by the rate d, and the kinetic rates of liga-
tion and cleavage are denoted by ki,j,L and kC respectively.
The arrays L and C represent the random graphs, chosen at
the beginning of each run. Using the symbol⊕ for sequence
concatenation, it holds thereby:
L(pi, pj , pk, pm) =
{
0 pi ⊕ pj *= pk
1with prob.rL pi ⊕ pj = pk
(3)
and
C(pi, pj , pk, pm) =
{
0 pi ⊕ pj *= pk
1with prob.rC pi ⊕ pj = pk
(4)
The sequence pm acts as catalyst.
In all calculations subsequently shown, several additional
assumptions have been made:
1. The monomers A,B must not act as catalysts; this in or-
der to enhance chemical plausibility.
2. There is a maximal sequence length L. Ligations may
well produce longer sequences, but those are assumed to
fall out by precipitation. This is physically plausible and
keeps the system tractable.
3. In order to capture steric effects, the ligation rate ki,j,L is
length dependent. Shall |pi| denote the length of pi, we
set ki,j,L = kL/(|pi||pj |) for some constant kL. The idea
behind this crude approximation is that in a well-stirred
reactor, the collision frequency of two sequences doesn’t
depend on the length. The collision happens by the con-
tact of two monomers, one out of each sequence. The
chance that those are the ones that are able of mutual liga-
tion because they mark the end and the start of the respec-
tive sequences is inversely proportional to the respective
length of the sequences.
The system then contains 2(L+1)−2 variables. This means,
taking into account the non-catalycity of the monomers, that
there are (2(L + 1) − 2)2(2(L + 1) − 4) potential ligation
reactions and (2(L+1)− 4)∑Ll=2 2l(l− 1) possible cleav-
age processes. As it turned out, already values of L = 6
deliver systems of sufficient combinatorial variety in order
to exhibit interesting dynamical effects. In all simulations,
we set ∀i : pi(0) = 1 as initial condition; this with the idea
to give a potential ACS in a random graph sufficiently favor-
able starting conditions. Following (13), a random reaction
graph qualifies as containing an ACS, if the concentration
of at least one non-monomeric species is above a threshold
T after more than 10td with td = − log(T )/d denoting the
typical decay time for T . As will be shown (and has already
be discussed by Bagley and Farmer), the decision whether a
reaction system contains an ACS as a subsystem is surpris-
ingly insensitive to the choice of T . The numerical solutions
were obtained by internal routines of the software package
MathematicaTM and a sample of solutions was verified with
a standard adaptive fourth-order Runge-Kutta solver.

!4 !2 2 4
log!kE "
5
10
15
20
25
30
35
average size ACS
Figure 7: Average size of the observed ACS in a random
reaction system with L = 6, kL = kC = 0.01 as a function
of the rate of energy influx kE and for a detection threshold
T = 10−2.
the situation of a flow reactor, there must not only be a path-
way for the production of a given molecule but its production
has in addition to happen at a rate that compensates for the
loss by outflow. Studying the kinetic behavior of random re-
action systems reveals the importance of a proper balancing
of the probabilities for different types of reaction: We inves-
tigated cleavage and found that taking into account dynam-
ics, cleavage does not only enlarge the variety of polymer
species (which is desirable from the perspective of obtain-
ing ACS) but may also destroy components relevant for the
system with a rate that cannot be compensated by their re-
spective generation processes. We also investigated the role
of energy consumption and found that the introduction of
energy as a limiting factor strongly influences the concen-
tration profile of the ACS. It turned out that whereas a large
supply of energy leads to a broad variability of sequences,
intermediate values seem to favor ACS with less, but, with
respect to concentration also in absolute terms, more pro-
nounced components.
We investigated systems with rather short sequences,
mostly with a maximal sequence length of L = 6. The nu-
merical values for the catalytic probabilities rL and rC need
then to be of a size which is chemically not realistic. We
claim, however, our results to be of worth because whereas
the quantitative features of the shown results heavily de-
pend on L, the qualitative don’t. Even more, data (partially
not shown) suggests that the discussed effects become more
pronounced with increasing L. According investigations
need then to be performed in a particle based manner, see
(3). Another interesting perspective is presently investigated
by DeLucrezia ((DeLucrezia and coworkers)) and cowork-
ers. In their approach, the “monomers” are replaced by pre-
prepared strands consisting of some ten amino acids. A se-
quence consisting of six of these strands may have a higher
probability of exhibiting catalytic properties. However, the
model presented in this paper is then only a “coarse-grained”
approximation to the dynamics, because cleavage may well
happen within one of the original monomeric strands.
Our choice of the initial conditions, namely to set the con-
centrations of all sequences under consideration to one at
the start is certainly unrealistic and motivated by our focus
on stability considerations. The discovery that the energy
supply influences the concentration profile opens the per-
spective of “iterative” emergence. A very limited set of ini-
tially provided components may establish a first, still frail
ACS which produces as side products some further, possi-
bly catalytic components at low concentrations. A temporal
increase of the energy supply may enable the system to reach
a new basin of attraction by a short-term increase of cleav-
ing catalysts which in turn produce a passing wider variety
of sequences. We will address this scenario in a subsequent
work focussed on issues of emergence, also considering as-
pects of stabilization against molecular parasites achieved
by spatial organization with (2) or without (5; 4) explicit
compartmentalization.
Taking into account dynamics shows that first, one of the
reasons for the fact that spontaneously formed autocatalytic
system have not or only rarely been observed in the labora-
tory may not only be due to lack of catalytic activity. As a
matter of fact, it could even be caused by too much cataly-
sis, if cleavage is too frequent. Second, and probably more
important, we need to shift our attention from focussing
solely on catalysis (and respective probabilities) to a picture
in which kinetics plays an important role too. Even if we had
reaction system in which in principle an ACS could produce
measurable signals, it only does if the kinetic parameters are
suitably chosen.
Acknowledgments. This work has been supported by the
Fondazione di Venezia, http://www.fondazionedivenezia.it,
(DICE project).
References
DeLucrezia, D. and coworkers. unpublished work.
Farmer, J. D., Kauffman, S. A., and Packard, N. H. (1986). Auto-
catalytic replication of polymers. Physica D, 22(2):50–67.
Filisetti, A., Serra, R., Carletti, T., Villani, M., and Poli, I. (2008).
Synchronization phenomena in protocell models. Biophysical
Reviews and Letters, 3(1/2):325–342.
Filisetti, A., Serra, R., Villani, M., Fu¨schslin, R. M., Packard, N.,
Kauffman, S. A., and Poli, I. (2010). A stochastic model of
autocatalytic reaction networks. submitted to ECCS10.
Fchslin, R. M., Altmeyer, S., and McCaskill, J. (2004). Evolu-
tionary stabilization of generous replicases by complex for-
mation. Europ. Phys J. B, 38(1):103–110.
Fchslin, R. M. and McCaskill, J. S. (2001). Evolutionary self-
organization of cell-free genetic coding. Proceedings of the
National Academy of Science USA, 98(16):9185–9190.

Hanel, R., Kauffman, S. A., and Thurner, S. (2005). Phase tran-
sition in random catalytic networks. Physical Review E,
72:036117.
Hordijk, W. and Steel, M. (2004). Detecting autocatalytic, self-
sustaining sets in chemical reaction systems. Journal of The-
oretical Biology, 227:451–461.
Jain, S. and Krishna, S. (2001). A model for the emergence of
cooperation, interdependence, and structure in evolving net-
works. Proceedings of the National Academy of Science USA,
98(2):543–547.
Kauffman, S. (1986). Autocatalytic sets of proteins. Journal of
Theoretical Biology, 119:1–24.
Letelier, J.-C., Soto-Andrade, J., Abarzua, F. G., Cornish-Bowden,
A., and Cardenas, M. L. (2006). Organizational invariance
and metabolic closure: Analysis in terms of (m;r) systems.
Journal of Theoretical Biology, 238:949–961.
Lifson, S. (1997). On the crucial stages in the origin of animate
matter. J. Molecular Evolution, 44:1–8.
Mossel, E. and Steel, M. (2005). Random biochemical networks:
the probability of self-sustaining autocatalysis. Journal of
Theoretical Biology, 233:327–336.
R.J.Bagley and J.D.Farmer (1991). Spontaneous emergence of a
metabolism. Artifcial Life II, pages 93–140. Santa Fe Institute
Studies in the Sciences of Complexity X.
Stadler, P. F., Fontana, W., and Miller, J. H. (1993). Physica D,
63:378.
Steel, M. (2000). The emergence of a self-catalysing structure in
abstract origin-of-life models. Appl. Math. Lett., 3:91–95.