Eur. Phys. J. B 77, 249–256 (2010)
DOI: 10.1140/epjb/e2010-00175-5
Regular Article
THE EUROPEAN
PHYSICAL JOURNAL B
Non-linear protocell models: synchronization and chaos!
A. Filisetti1,a, R. Serra1,2, T. Carletti1,3, M. Villani1,2, and I. Poli1,4
1 European Centre for Living Technology, Calle del Clero 2940, 30124 Venezia, Italy
2 Dipartimento di Scienze Sociali, Cognitive e Quantitative, Universita` di Modena e Reggio Emilia, via Allegri 9,
42100 Reggio Emilia, Italy
3 De´partement de mathe´matique, Faculte´s Universitaires Notre Dame de la Paix, rempart de la Vierge 8,
B 5000 Namur, Belgium
4 Dipartimento di Statistica, Universita` Ca’ Foscari, Giobbe - Cannaregio 873, 30121 Venezia, Italy
Received 13 December 2009 / Received in final form 22 April 2010
Published online 10 June 2010 – c© EDP Sciences, Societa` Italiana di Fisica, Springer-Verlag 2010
Abstract. We consider generic protocells models allowing linear and non-linear kinetics for the main in-
volved chemical reactions. We are interested in understanding if and how the protocell division and the
metabolism do synchronise to give rise to sustainable evolution of the protocell.
1 Introduction
Protocells could be lipid vesicles, or micelles, endowed
with some rudimentary metabolism and should contain
“genetic”material, being able to grow, reproduce and
evolve. While viable protocells do not yet exist, their study
is important in order to understand possible scenarios
for the origin of life, as well as for creating new “pro-
tolife” forms which are able to adapt and evolve [1]. This
endeavour has an obvious theoretical interest, but it might
also lead to an entirely new “living technology”, definitely
different from conventional biotechnology.
Theoretical models can be extremely useful to devise
possible protocell architectures and to forecast their be-
haviour, a review of the current knowledge is available
in [2,3]. What can be called the “genetic material” of a
protocell is composed by a set of molecules which, col-
lectively, are able to replicate themselves. At the same
time, the whole protocell undergoes a growth process (its
metabolism) followed by a breakup into two daughter
cells. This breakup is a physical phenomenon which is
frequently observed in lipid vesicles, and it has nothing
to do with life, although it superficially resembles the di-
vision of a cell [4]. In order for evolution to be possible,
some genetic molecules should affect the rate of duplica-
tion of the whole container, and some mechanisms have
been proposed whereby this can be achieved.
In this paper we address an important issue in proto-
cell research: the genetic material duplicates at a certain
rate, while the lipid container grows, in general, at an-
other rate. When the container splits into two, it may
! This work has been originally presented at the European
Conference on Complex Systems, held in Warwick from 21 to
25 September 2009.
a e-mail: alessandro.filisetti@ecltech.org
happen that the genetic material has not yet doubled: in
this case its concentration would be lower in the daughter
protocells. Hence through generations, this density might
eventually vanish. On the other hand, if the genetic ma-
terial growth were faster than the container, the former
would accumulate in successive generations.
So, in order for a viable population of evolving pro-
tocells to form, it is necessary that the rhythms of the
two processes are synchronised. In some models (like the
Chemoton [5] this is imposed a priori in the kinetic equa-
tions, but it is unlikely that such a set of exactly coupled
reactions springs up spontaneously in a single step. It is
therefore interesting to consider the possibility that such
synchronization could be an emergent phenomenon, with-
out imposing it a priori.
It has previously been shown that this may indeed hap-
pen when one takes into account successive generations of
protocells [6]. Even if at the beginning the two processes
take place with different paces, they asymptotically ap-
proach synchronization, in the sense that both (i) the time
interval between the formation of a protocell and the mo-
ment when it doubles its size and (ii) the time required for
doubling the genetic material change, do vary and gener-
ation after generation, they tend to the same value.
The present paper explores under which conditions
such synchronization takes place, in the interesting case
where there are several different kinds of replicating
molecules in each protocell, and they can affect each
other’s replication rates in a non-linear way. It therefore
expands our previous studies, which had considered syn-
chronization (i) when there is only one kind of replicating
molecule in each protocell [6] (ii) when there are differ-
ent kinds of replicators in each protocell, but they do not
directly affect each other’s replication rate [6] and (iii)
when there are different kinds of replicators in each pro-
tocell, which affect each other’s replication rates, and the

250 The European Physical Journal B
Fig. 1. A schematic representation of the Surface Reaction
Models. In the figure a portion of the vesicle membrane is
shown.
kinetic equations are linear [7]. Note that the model as
a whole is always non-linear, due to the division of the
lipid container: linearity refers only to the set of kinetic
equations for the replicators in a protocell.
Synchronization is studied here using abstract mod-
els belonging to the “surface reaction models” family
(briefly, SRMs); let us observe that other architectures
have been introduced, for instance the internal reaction
models (IRMs) [7]. The difference is that, in the former
case, the reactions which lead to the formation of the new
genetic material and those devoted to the formation of
new membrane molecules take place close to the proto-
cell outer surface, while in IRMs they both take place in
the interior of the vesicle. The modelling level is fairly
abstract, so the results should hold for different detailed
protocell architectures. SRMs are inspired by the the so-
called “Los Alamos bug”, a model of protocells where the
genetic material is composed by strands of peptide nu-
cleic acids, PNA for short [8,9], which should be found in
the vesicle membrane. According to this hypothesis, dif-
ferent PNA’s may influence the growth rate of their “con-
tainer” by catalysing the formation of amphiphiles (which
form the protocell membrane) from precursors, Figure 1.
On the other hand IRMs are supported among other, by
the RNA-cell architecture proposed by [10–12].
The paper is organised as follows. For the sake of com-
pleteness in Sections 2 and 3 we recall the main features
of SRMs, the main mathematical techniques to study syn-
chronization and the results of our previous studies. Sec-
tions 4, 5 and 6 are devoted to the introduction of some
non-linear replicators kinetic models and the description
of their synchronization properties. Section 7 is aiming at
understanding the outcome of regular behaviours once the
involved kinetic equations exhibit chaotic regimes: embed
chaos to make it disappear. Finally, in the last section
some critical comments and indications for further work
are reported.
2 Surface reactions models
Let us suppose the protocell to be a spherical vesicle, with
an aqueous interior and a lipid-phase membrane composed
by amphiphilic molecules. In a way inspired by the Los
Alamos bug hypothesis, we will suppose that the mem-
brane grows by addition of amphiphiles, which are formed
close to its external surface by suitable precursors, under
the catalytic influence of some molecules which are found
in the vesicle, hereby called Genetic Memory Molecules,
GMMs for short. In general, only a fraction of these cata-
lysts will be effective, namely those which are close enough
to the outer surface (it is assumed that these molecules can
be found in the lipid phase). We will also suppose that the
membrane is composed by a single kind of amphiphilic
molecules, Figure 1.
Let us then consider a single protocell, and let C be the
quantity of its lipid constituent (moles of amphiphiles).
Then the growth of C in time is assumed to be propor-
tional to the vesicle surface S, times a function of the
concentration of the catalysts in the membrane, times a
function of the density of the precursors which are found
outside, close to the membrane. We will assume that this
last term is not influential, e.g. that either precursors are
buffered, so that their concentration is kept constant, or
that there are saturation effects and the concentration of
the precursors is high enough to saturate.
Moreover, we will also assume (i) that the rate of spon-
taneous formation of amphiphiles in the outer medium is
negligible with respect to the one due to the catalytic ef-
fects of the replicators and (ii) that diffusion is fast enough
to ensure that, on the time scale of the model, the cata-
lyst concentration is homogeneous in the lipid phase: so,
if there is a single kind of catalyst, and if X denotes its
quantity (moles) in the protocell lipid phase, its concen-
tration is [X ] = X/VL, where VL is the volume of the lipid
phase.
Therefore, under the above hypotheses
dC
dt = S(VL)fˆ
(
X
VL
)
. (1)
Where fˆ denotes a function of the concentration that will
be specified later on, and the dependency of S upon the
volume VL has been emphasised. We will assume that S is
proportional to V βL , with the parameter β ranging in the
interval [2/3, 1]: in a spherical micelle, where the volume
of the lipid phase equals the total volume, S would be
proportional to V 2/3L , while in a vesicle with a very thin
membrane it would be almost proportional to VL itself.
Since VL = C/ρ (with ρ constant), by redefining the arbi-
trary function which appears in the r.h.s. of the previous
equation we obtain
dC
dt = C
βf
(
X
C
)
. (2)
If, as it is often the case, we assume f to be proportional
to X/C to some power γ we obtain
dC
dt = αC
β−γXγ . (3)
A particularly important studied case is when the growth
of the container is assumed to be proportional to the

A. Filisetti et al.: Non-linear protocell models: synchronization and chaos 251
concentration of the catalyst molecules, i.e.
dC
dt = αC
β−1X. (4)
Equations (3) and (4) describe the continuous growth of
the protocell up to a certain size, which is then followed
by the breakup in two daughter cells. Although this is a
complicated process, we will assume for simplicity (as it is
done in most models) that the cell breaks when it reaches
a certain size, i.e. when C(t) equals a fixed value θ and
that the initial value of C of each of the two daughter cells
is θ/2.
We also suppose, like in Los Alamos bug, that also the
replication of the X -type molecules takes place near the
external surface, by virtue of buffered precursors which are
also in this case not taken explicitly into account. Replica-
tion takes place in a thin boundary of width δ near the sur-
face, therefore only the molecules which are found in this
volume can contribute to the synthesis of new molecules.
Therefore dX /dt should be proportional to Sδ times a
function of X/C , so
dX
dt = S(V
β
L )gˆ
(
X
C
)
= Cβg
(
X
C
)
. (5)
Assuming g to be proportional to X/C raised to some
power ν we then obtain
dX
dt = ηC
β−νXν. (6)
The case where the continuous growth phase is described
by equations (4) and (6) will be referred to as the “single
replicator”case.
Note that in deriving equations (5) and (6) it has been
implicitly assumed that replication needs the presence of
externally supplied precursors, which are available only
close to the outer boundary of the membrane. This is
consistent with a template duplication mechanism, like
that of nucleic acids; if replicators were e.g. interacting
polypeptides, and precursors were available in the whole
lipid phase, then all the X molecules in the lipid phase
could contribute and, instead of equation (5) we would
get
dX
dt = Cg
(
X
C
)
. (7)
And, if g is proportional to X/C raised to some power ν
dX
dt = ηC
1−νXν . (8)
So the difference between the two cases would amount
to choosing β = 1 in the second one. But it has already
been shown [6] that, while the kinetics is influenced by
the value of the geometric parameter β, the achievement
of synchronization is not, so that one can limit to analyse
the case β = 1 to draw general conclusions. We therefore
come to the conclusion that, as far as synchronization is
concerned, both a “nucleic-acid based” and a “polypeptide
based” hypothesis are described by the formalism used
here and therefore behave in much the same way.
A particularly important case, whose analysis is quite
simple, is the one where both exponents γ and ν are equal
to one, so {
dC
dt = αCβ−1X
dX
dt = ηCβ−1X.
(9)
Which can be straightforwardly generalised to the case of
N different linearly interacting replicators :
{
dC
dt = αXCβ−1
dX
dt = MXCβ−1.
(10)
Where α and X are N-dimensional vectors and M is the
N ×N matrix of kinetic coefficients. The same arguments
used above concerning the role of β, we can limit ourselves
to consider the case β = 1 in the study of the synchro-
nization phenomenon.
Next step is to consider two replicators, X and Y , in
the same protocell whose growth rate is assumed to be
non-linear. If they do not interact directly, and their effect
on the growth of the lipid container is linear, then the
continuous growth is described by the following equations,
analogous to equation (6)



dC
dt = α′X + α′′Y
dX
dt = η′XνCβ−ν
dY
dt = η′′Y νCβ−ν .
(11)
A more interesting case is when the replicators directly
interact each others. For instance, let us suppose that the
kinetics is second order, so that the rates of production
of new X and new Y should be proportional to the fre-
quency of encounters between X and Y molecules. The
total number of encounters per second, in the lipid phase,
is then Sδ[X ][Y ] ≈ V β−2L δXY . By observing that VL is
proportional to C and by lumping some constants in the
terms η′ and η′′ one then gets



dC
dt = α′Cβ−1X
dX
dt = η′Cβ−2XY
dY
dt = η′′Cβ−2XY.
(12)
Where we have assumed for simplicity that only one kind
of molecule, say X , catalyses the growth of the membrane,
i.e. set α′′ = 0 in equation (11). In the more general case
where the rate of production of X is proportional to the
concentration of Y times some power ν ofX , and similarly
for Y , the previous equations generalise to



dC
dt = α′Cβ−1X
dX
dt = η′Cβ−1−νXνY
dY
dt = η′′Cβ−1−νXY ν .
(13)
Equations (11) and (12) are examples of cases with two
non-linear interacting replicators. More general non-linear
equations will be introduced in Section 4. Let us remark
that even if the quadratic model (12) is widely used in
literature, see for instance [13,14], it doesn’t generically
achieve the synchronization.

252 The European Physical Journal B
3 A summary of previous results
The equations of SRMs lend themselves to a nice analyt-
ical technique, which has been applied in order to study
their behaviour, and which has been complemented by nu-
merical simulations1.
Starting with an initial quantity of container C at time
T0 equal to θ/2, we assume that once C reaches a critical
value θ, the protocell will divide into two equal protocells
each one of size θ/2. Let ∆T0 be the time interval needed
to double C from this initial condition, and let T1 = T0 +
∆T0 be the time when the critical value θ is reached. Since
the initial value for C is fixed, ∆T0 is a function of the
initial quantity of replicators,X0. The final value ofX in
the protocell, just before the division will be denoted by
X(T1) . Because we assume perfect halving at the division,
each offspring will start with an initial concentration of
replicators equal to X1 =X(T1)/2. Then the continuous
growth process starts again, until the next doubling time
will be reached, T2 = T1 + ∆T1, and the third generation
will start with an initial value X2 =X(T2)/2, and so on.
We generalise the preceding discussion with the fol-
lowing equations, which refer to the kth cell division cycle
that starts at time Tk and ends at time Tk+1:
θ
2 =
∫ Tk+1
Tk
dC
dt dt , and Xk+1 =
1
2X(Tk+1). (14)
Note that in general X(Tk+1) #= 2X(Tk) and therefore
the time needed to double the value of C is not constant
between two successive generations. The problem of syn-
chronization therefore amounts to find under which con-
ditions the following equality asymptotically holds:
lim
t→∞
[X(Tk+1)− 2X(Tk)] = 0. (15)
The existence of first integrals of the equations describ-
ing the continuous growth phase, allows us to prove syn-
chronization. In fact we are able to obtain a discrete
map X(Tk+1) = F (X(Tk)), and the search for condi-
tions which allow synchronization are equivalent to the
existence of a fixed point.
There is no general form for the first integrals, which
depend upon the kind of SRM model which is consid-
ered: several examples can be found in [6,7,15] and the
main results are hereby summarised. Before doing so, it is
however interesting to consider which kind of behaviours
one can imagine to find. The possible alternatives are the
following:
1 Synchronization: in successive generations (as k → ∞)
the time for duplication of the protocell, ∆Tk, and the
time required to duplicate the genetic material, ∆T gk ,
approach the same value (for the reasons highlighted
1 Numerical simulations have been performed using Matlab’s
standard solver for ordinary differential equations ode45 with
parameter “nonNegative”. This function implements a Runge-
Kutta method with a variable time step for efficient computa-
tion and prevents the variables to become negative, for details
see Matlab website (http://www.mathworks.com).
above, this condition is satisfied when equation (15)
holds).
2 The initial concentration of the genetic material van-
ishes in the limit of infinitely many divisions; in this
case, given the above assumptions, the growth of the
container ends and the whole process stops.
3 The initial concentration of the genetic material grows
unbounded a k → ∞ . This results in a limitation of
the previous equations to modelling a protocell, which
indeed lack a rate limiting term for the growth rate of
X .
4 The two rates ∆Tk and ∆T gk oscillate in time with the
same frequency. This condition is not equivalent to
synchronization strictu sensu but it would nonethe-
less allow sustainable growth of the population of
protocells. Therefore this condition might be called
super-synchronization. Note that in principle super-
synchronization does not require equality of the two
frequencies, but that their ratio be a rational number.
6 The two rates ∆Tk and ∆T gk exhibit a non regular time
behaviour.
In the case of a single replicator whose replication rate
given by equation (6) with 0 < ν < 2, with linear growth
for the container, one can show that synchronization is
always achieved provided that α and η are positive (these
conditions are fairly obvious, otherwise it would make no
sense to speak of growth). The result holds also if one
takes a realistic geometry into account (i.e. if one uses the
volume of a spherical shell instead of the approximation
S ≈ V βL ). It also generalises to more general forms of the
function which describes the growth of the protocell or
that of the replicators [16]. More generically, we can prove
that once the container growth follows some non-linear
power law, i.e. γ #= 1 in equation (4), then a sufficient
condition to have synchronization is 0 < ν < γ + 1.
Let us observe that the interesting case of parabolic
growth, i.e. ν < 1, fits the above condition. On the other
hand, we can prove and it has been numerically checked
that in the case of a very steep increase ν ≥ 2, the system
does not approach synchronization: if ν = 2 then synchro-
nization can be obtained only with a very specific choice of
the involved parameters, whereas if ν > 2 synchronization
can never be achieved, according to the initial amount of
X , in the long run either Xk diverges or goes to zero.
A further generalisation is worth discussing. In
Section 2 we have assumed that [X ] = X/VL, which is
correct as long as X itself does not appreciably contribute
to the volume of the lipid phase. But if the quantity of
X becomes large, and X itself is a lipophilic compound
which contributes to the container, this formula should be
substituted by [X ] = X/(VL+VX). This can be rewritten,
by rescaling the kinetic constants, as X/(C + rX). It has
been analytically shown and numerically confirmed that
also in this case synchronization is achieved when α and
η are positive.
If there are several replicators in the same cell, but
they do not interact directly, namely the system is mod-
elled by equation (11), one again finds synchronization if
α, η and η′ are positive. In the linear case, ν = 1 , with two

A. Filisetti et al.: Non-linear protocell models: synchronization and chaos 253
replicators, only the fastest replicator survives in the final
population of protocells, while if ν < 1 both survive, their
relative proportion being a function of η/η′ . This is con-
sistent with similar behaviours observed in population dy-
namics. On the other hand if there is mutual interaction,
namely the protocell can be described by equation (12)
or equation (13), then synchronization cannot be always
achieved.
The case of several linearly interacting replicators
equation (10) has been analytically studied and a com-
plete discussion of the different cases can be found in [7].
The most relevant results are that the behaviour of the
system is ruled, in the long time limit, by the eigenvalue
of the matrix M with largest real part, denoted by λLRP .
If )(λLRP) > 0, and if the corresponding eigenvector v1
is non-negative2, synchronization is achieved. The asymp-
totic value of Xk is a multiple of v1. A sufficient condi-
tion to guarantee that λLRP is real and positive and that
the associated eigenvector is non-negative, is that the ma-
trix elements Mij are non-negative. This is an important
case, corresponding to all X molecules, which directly in-
teract, contribute to the synthesis of the others (mutual
catalysis).
It is however possible to imagine also cases where the
network of reactions includes some negative terms (if they
were all non-positive the system would of course die out).
The most relevant phenomenon discovered in this analysis
is the existence of an oscillatory behaviour, supersynchro-
nisation, which is found when the eigenvalues with the
largest real part are complex conjugate (see for instance
Fig. 2).
Let us observe that the most striking result of the anal-
ysis of linear replicators is that they behave in a way simi-
lar to that of a Continuously Stirred Tank Reactor : in the
present case vesicle splitting limits the asymptotic values,
while in CSTR it is the outflow which does it. But the ra-
tio of various replicators is the same in the two cases. The
coefficient α does not influence this ratio, nor the asymp-
totic division time, although it affects the actual value of
the asymptotic quantities.
These latter cases show that non-linearity may lead
to a halting of the growth, or to an unbounded growth
of the molecules, hence preventing the protocell from a
viable evolution. These behaviours can be associated with
the presence of power laws growth rates, to tackle the
question we thus consider more general non-linear models
by introducing “squashing functions”.
4 Non-linear models with squashing terms
The growth of the container is still assumed to be de-
scribed by a linear equation as equation (4), for instance
dC
dt = αXC
β−1. (16)
2 Because eigenvectors are defined up to a multiplicative fac-
tor, by non-negative we mean that all the components can be
made positive.
0 20 40 60 80 1000,8
0,9
1
1,1
1,2
Division number (k)
∆
T
k
(a)
0 20 40 60 80 1000
50
100
150
200
Time at each division
X
k


(b)
Fig. 2. An example of 4 × 4 matrix M with negative en-
tries where synchronization in not achieved but the division
time and the initial amount of genetic material regularly os-
cillate division after division. Upper panel: the division time,
∆Tk, in function of the generation number k. Bottom panel:
the time evolution behaviour of the amounts of Xk at the be-
ginning of each generation. Symbols refer to X1 = #,X2 =
!,X3 = ◦,X4 = ♦ (C0 = 500, X0 = [20.3, 20.2, 27, 47], M
= [0, −0.0690, 1.6821, −3.7767; 1.6905, 0, −2.3801, 1.3082;
1.3082, 3.9692, 0, 4.283; 0.6255, 3.838, −4.1488, 0]).
We then select some possible kinetic equations for the
replicators growth, that we now briefly describe and anal-
yse with the help of dedicated numerical simulations.
4.1 Quasilinear models
One drawback of linear equations like equation (10) is that
the growth rate may undergo an unrealistic unlimited in-
crease. In order to take physical constraints on the rates
into account it is possible to introduce bounds which are
never exceeded. Instead of fixing sharp thresholds, which
would lead to discontinuous derivatives, we consider here
“squashing functions”, i.e. functions which are bounded
both from below and from above, and which are never de-
creasing. They can be though as activating functions for
the reactions to start.
Let σ(x) be such a function, then we hypothesise
the following replicator law in the case of N interacting
molecules:
dX i
dt = C
βσ
(
N∑
k=1
MikXk
C
)
. (17)
In this case the behaviour is similar to that of the cor-
responding linear model, i.e. equation (10), once we use

254 The European Physical Journal B
the same coefficients Mij . In particular, synchroniza-
tion in the linear models implies synchronization in the
quasi-linear one. However sometimes we have observed
super-synchronization in the linear case, while the cor-
responding non-linear version still synchronises. We fre-
quently observe that the amounts of X -molecules rapidly
saturate, so that the duplication times are largely unaf-
fected by the exact values of the matrix elements Mij and
are not a function of λLRP , as it would happen in the lin-
ear case. A second interesting observed feature is that cell
duplication times are not affected by the value of α.
5 Second order kinetics
The quadratic model equation (12) can be generalised as
to include N interacting GMMs, to give
dX i
dt = C
β−2
N∑
k=1
MikXiXk. (18)
In the case of mutual catalysis the matrix coefficients are
non-negative. Note that catalytic cycles can be modelled
in this way by a proper choice of the matrix elements Mik.
Generically this model does not show synchronization, ex-
cept for very peculiar relationships among the coefficients.
On the other hand, once we introduce molecules able
to self-replicate
dX i
dt = C
β−2
N∑
k=1
MikXiXk + Cβ−1ηiXi, (19)
we can have synchronization. Let us observe that fixing
non-negative matrix elements Mik in a cyclic way one gets
the well known hyper-cycle systems [17].
The fixed point corresponding to synchronization de-
pends on the value of η as numerically shown in Figure 3,
where we present two set of simulations: in the first one,
top panel, using a fixed value of η the system approaches
the same fixed point although it started from different
initial conditions; in the second one, bottom panel, four
different fixed points are reached using different values of
η, using the same initial conditions.
6 Second order kinetics without
self-replication
Let us finally consider the case where there is no self-
replication, but the GMMs mutually catalyse each other’s
formation from existing precursors, in a way it requires
the interaction of two molecules to produce a third one.
The corresponding equations, neglecting possible satura-
tion effects, are then
dXi
dt
= Cβ−2
N∑
k=1
MijkXjXk. (20)
0
500
1000
0
20
40
60
0
100
200
300
400
500
X1X2
X3
(a)
0
500
1000
0
500
1000
0
200
400
600
800
X1X2
X3
(b)
Fig. 3. Synchronization fixed point for the system ruled by
equation (19). On the top panel four different runs with fixed
values for η and different X0 are shown ( C0 = 500, X10 =
[2.77, 32.93, 17.72], X20 = [15.05, 33.49, 39.31], X
3
0 = [38.50,
42.48, 18.03], X40 = [10.60, 26.95, 19.06], η = [1.5, 0.5, 1] ). On
the bottom panel four different runs with random values of η
are shown( C0 = 500, X10 = [38.41, 7.29, 26.85], X
2
0 = [1.18,
29.93, 48.22], X30 = [13.45, 37.46, 16.29], X
4
0 = [21.29, 11.07,
15.81], η1 = [0.4139, 1.665, 1.214], η2 = [0.927, 0.432, 0.826],
η3 = [1.967, 1.167, 1.813], η4 = [1.276, 0.278, 0.034], ). In both
cases α = [1, 1, 1] and M = [0, 1, 0; 0, 0, 1; 1, 0, 0].
Where, in order to avoid self-replication, the matrix ele-
ments should have the form
Mijk = µijk(1− δij)(1− δik). (21)
In this case we observe that sometimes synchronization is
achieved, while in other cases extinction is the outcome.
This result is related to the sparseness of the matrix M . It
is interesting to observe that, if a large fraction of the ma-
trix elements is non-vanishing, synchronization is found,
even when the frequency of negative Mijk’s equals that of
positive ones.
Increasing the number of GMMs the relative number of
non zero entries increases, hence the matrix is less sparse.
Several simulations show that, considering an equal fre-
quency of negative and positive entries, ten replicators
(that correspond to a sparsity coefficient equal to 0.28)
are sufficient in order to avoid death of the protocell.

A. Filisetti et al.: Non-linear protocell models: synchronization and chaos 255
0
50
100
0
50
100
0
10
20
30


Fig. 4. The Willamowsky-Ro¨ssler strange attractor (X(0) =
2, Y(0) = 3, Z(0) = 1, k1 = 30, k2 = 1, k3 = 10, k4 = 1, k5 =
16.5, k−1 = 0.25, k−2 = 0.0001, k−3 = 0.001, k−4 = 0.5,
k−5 = 0.5).
7 Using a chaotic kinetic for the genetic
memory molecules
The models considered so far are, from a mathemati-
cal point of view, generic non-linear systems of differ-
ential equations, that we know can exhibit chaotic be-
haviours. It is thus a natural question to consider how
these non-regular behaviours can be associated to sustain-
able growth and evolution of protocells. In other words,
why, even if you have a set of chemical reactions that
exhibit a chaotic behaviour, once you “embed” these re-
actions in a protocell, their behaviour is more regular, in
such a way the protocell can grow and divide?
In the following we provide a partial answer to this
question by considering the Willamowsky-Ro¨ssler [18] sys-
tem3, that exhibits chaotic behaviours (see Fig. 4) and it
has been already used to model chemical reactions:



dX
dt = k1X − k−1X2 − k2XY + k−2Y 2
−k4XZ + k−4
dY
dt = k2XY − k−2Y 2 − k3Y + k−3
dZ
dt = −k4XZ + k−4 + k5Z − k−5Z2.
(22)
We observe that if the previous kinetic equations are intro-
duced in a protocell, whose container varies as a function
of some GMMs, then the Willamowsky-Ro¨ssler protocell
becomes, under the assumption β = 1:



dC
dt = αX
dX
dt = k1X − 1C k−1X2 − 1C k2XY
+ 1C k−2Y 2 − 1C k4XZ + Ck−4
dY
dt = 1C k2XY − 1C k−2Y 2 − k3Y + Ck−3
dZ
dt = − 1C k4XZ + Ck−4 + k5Z − 1C k−5Z2.
(23)
We observed that the system undergoes through a bi-
furcation as the parameters α varies: for small values of
3 Similar conclusions have been found by analysing a mod-
ified Lorentz system where the strange attractor is fully con-
tained in the positive octant.
0
1000
2000
3000
0
500
1000
1500
0
5
10
15

X1X2

X
3
(a)
0
2000
4000
0
500
1000
1500
0
200
400
600

X1X2

X
3
(b)
0
2000
4000
0
2000
4000
0
500
1000
1500

X1X2

3
(c)
Fig. 5. Graphs show the evolution of the system changing the
value of parameter α. For a small value of α (a) the chaotic
behaviour is maintained but increasing the latter the system
reaches a periodic orbit at first (b) and then synchronization
(c) (Parameters: α: 0.01(a), 0.1(b) and 0.3(c), X(0) = 2, Y(0) =
3, Z(0) = 1, k1 = 30, k2 = 1, k3 = 10, k4 = 1, k5 = 16.5, k 1 =
0.25, k 2 = 0.0001, k 3 = 0.001, k 4 = 0.5, k 5 = 0.5).
the latter the system still evolves chaotically following
the Willamowsky-Ro¨ssler strange attractor, then as α in-
creases the system presents, first a periodic orbit and then
a fixed point, hence synchronization, see Figures 5 and 6.
This implies that protocells, and possibly also real cells,
during their evolution could have “tuned”, or have been
selected, in such a way the coupling between container
growth and metabolism/information reproduction, gives

256 The European Physical Journal B
0,1 0,09 0,08 0,07 0,06 0,05 0,04 0,03 0,02 0,01

0
1
2
3
4
5
6
7
T
c
Fig. 6. The behaviour of the Willamowsky-Ro¨ssler system in-
troduced in a protocell is tightly correlated with the value of α.
The figure on the top panel represents the bifurcation diagram
of ∆TC in function of α variation. On the X axis α values are
represented while on the Y axis the cell division time is shown
(α goes from 0.1 to 0.01 and for each value of the latter the last
twenty cell division times are shown). Although it is not clear
observing the graph, the system shows supersynchronisation
for α equal to 0.1. On the bottom panel Z in function of X
is represented using different values for the parameter α. It is
clearly observable that the orbits dimension decrease increas-
ing the value of α approaching the synchronization fixed point.
For values of α larger than 0.3 (here not shown) the limit cy-
cle became a fixed point and synchronization is achieved (the
values of the parameters are the same of Fig. 4).
rise to a regular behaviour, even if the latter systems could
evolve (separately) in a chaotic way.
8 Conclusions
In this work studies on linear reaction models have been
extended to the non linear ones. Excepting for the second
order kinetic model, Section 5, the emergent synchroniza-
tion phenomena are maintained showing a high degree
of robustness of these models against perturbations in
kinetics. Interesting outcomes have been obtained cou-
pling the container kinetic with a chaotic set of interacting
molecules showing that the evolution of the system is re-
lated to the growth rate of the protocell: the more the
growth is fast, the more the protocell is able to grow and
divide. The addressed questions are surely relevant to un-
derstand the emergence of sustainable forms of life, and
thus deserves further investigations that will be presented
elsewhere.
Authors would like to thank the “DICE project of Fondazione
Venezia”and the European Centre for Living Technology for
the financial and logistic support.
References
1. S. Rasmussen, L. Chen, D. Deamer, D.C. Krakauer, N.H.
Packard, P.F. Stadler, M.A. Bedau, Science 303, 963
(2004)
2. R.V. Sole´, J. Mac´ıa, H. Fellermann, A. Munteanu, J.
Sardanye´s, S. Valverde, 213 (2008)
3. R.V. Sole´, A. Munteanu, C. Rodriguez-Caso, J. Mac´ıa,
Synthetic protocell biology: from reproduction to compu-
tation. Philos. Trans. R. Soc. Lond. B Biol. Sci. 362, 1727
(2007)
4. M.M. Hanczyc, J.W. Szostak, Curr. Opin. Chem. Biol. 8,
660 (2004)
5. T. Ganti, Chemoton theory: Theory of fluid ma-
chineries: Theory of living system (New York, Kluwer
Academic/Plenum, 2003), Vol. i, ii
6. R. Serra, T. Carletti, I. Poli, Artificial Life 13, 1 (2007)
7. T. Carletti, R. Serra, I. Poli, M. Villani, A. Filisetti, J.
Theor. Biol. 254, 741 (2008)
8. S. Rasmussen, L. Chen, M. Nilsson, S. Abe, Artificial Life
9, 269 (2003)
9. S.J. Rasmussen, L. Chen, B.M.R. Stadler, P.F. Stadler,
Orig. Life Evol. Biosph. 34, 171 (2004)
10. S.S. Mansy, J.P. Schrum, M. Krishnamurthy, S. Tobe´, D.A.
Treco, J.W. Szostak, Nature 454, 122 (2008)
11. T. Oberholzer, R. Wick, P.L. Luisi, C.K. Biebricher,
Biochem. Biophys. Res. Commun. 207, 250 (1995)
12. J.W. Szostak, D.P. Bartel, P.L. Luisi, Nature 409, 387
(2001)
13. K. Kaneko, Adv. Chem. Phys. B 130, 543 (2003)
14. K. Kaneko, Life: An Introduction to Complex Systems
Biology (Understanding Complex Systems) (Springer-
Verlag New York, Inc., Secaucus, NJ, USA, 2006)
15. A. Filisetti, R. Serra, T. Carletti, M. Villani, I. Poli,
Biophys. Rev. Lett. (BRL) 3, 325 (2008)
16. R. Serra, T. Carletti, I. Poli, M. Villani, A. Filisetti,
Conditions for emergent synchronization in protocell, in
Proceedings of ECCS07: European Conference on Complex
Systems, edited by J. Jost, D. Helbing CD-Rom, paper No.
68, 2007
17. M. Eigen, P. Schuster, Naturwissenschaften 64, 541 (1977)
18. B.D. Aguda, B.L. Clarke, J. Chem. Phys. 89, 7428 (1988)