hr
ise
art d
1 Ve
gio E
Received 22 January 2008
Received in revised form
Available online 16 July 2008
Keywords:
Protocell
Self-replication
Dynamical model
Synchronization
g
phenomenon of genetic material and container productions, a necessary condition to ensure sustainable
growth in protocells and eventually leading to Darwinian evolution when applied to a population of
able to capture universal behaviors, without carefully adding
or amphiphile) and one kind of replicator molecule—loosely
rates,
tocell
n, to
than duplication, the concentration of genetic material would
ARTICLE IN PRESS
Contents lists availabl
.e
Journal of Theor
Journal of Theoretical Biology 254 (2008) 741–751genetic material at duplication time refers to the average beha-Tel.: +32 81724903; fax: +3281724914.E-mail addresses: timoteo.carletti@fundp.ac.be (T. Carletti), rserra@unimore.iteventually vanish (we refer to the splitting of a protocell as
duplication, and to the doubling of genetic polymers as replica-
tion). In the opposite case, its concentration would grow
unbounded. Of course, the requirement of doubling of the
vior, while each single event is affected by noise and random
fluctuations.
 Corresponding author at: De´partement de Mathe´matique, Faculte´s Universi-
taires Notre Dame de la Paix, Rempart de la Vierge 8, B 5000 Namur, Belgium.
(R. Serra), irenpoli@unive.it (I. Poli), villani.marco@unimore.it (M. Villani),
alessandro.filisetti@ecltech.org (A. Filisetti).0022-51
doi:10.1complicating details’’ (Kaneko, 2006). A protocell should comprise
at least one kind of ‘‘container’’ molecule (typically a lipid
ensure that each offspring will contain the same amount of
genetic material as the mother. Indeed, if replication were slowerIn order to study how protocells can develop, given that they
do not yet exist, it is necessary to consider ‘‘simplified models
by dilution, it is necessary that the two proceed at equal
i.e., that the genetic material has doubled when the pro
splits into two—a condition referred to as synchronizatioSeveral attempts are currently under way to obtain protocells
capable of growth and duplication, endowed with some limited
form of genetics (Oberholzer et al., 1995; Rasmussen et al., 2004;
Szostak et al., 2001; Mansy et al., 2008). The interest for these
systems is motivated either by the quest to understand which are
the minimal requirements for a life form to exist and evolve, or by
the search for indications about the way in which primitive life
might have emerged on Earth.
can be copied or a system of two or more kinds of replicators
which catalyze each other’s synthesis—e.g., proteins and nucleic
acids. There are therefore two kinds of reactions which are crucial
for the working of the protocell, which in this paper will be called
‘‘key’’ reactions: those which synthesize the container molecules
and those which synthesize the GMM replicators.
The two key reactions may take place at different rates.
However, to achieve sustained protocell growth and avoid death1. Introduction93/$ - see front matter & 2008 Elsevier Ltd. A
016/j.jtbi.2008.07.008Synchronization has been proved to be an emergent property in many relevant protocell models in
the class of the so-called surface reaction models, assuming both linear- and non-linear dynamics for
the involved chemical reactions. We here extend this analysis by introducing and studying a new class
of models where the relevant chemical reactions are assumed to occur inside the protocell, in contrast
with the former model where the reaction site was the external surface.
While in our previous studies, the replicators were assumed to compete for resources,
without any direct interaction among them, we here improve both models by allowing linear
interaction between replicators: catalysis and/or inhibition. Extending some techniques previously
introduced, we are able to give a quite general analytical answer about the synchronization
phenomenon in this more general context. We also report on results of numerical simulations to
support the theory, where applicable, and allow the investigation of cases which are not amenable to
analytical calculations.
& 2008 Elsevier Ltd. All rights reserved.
speaking ‘‘genetic material’’, hereafter called, genetic memory
molecule, GMM for short. This is typically a linear polymer which8 July 2008
Accepted 8 July 2008
protocells.Sufficient conditions for emergent sync
T. Carletti a,b,, R. Serra b, I. Poli b, M. Villani c, A. Fil
a De´partement de Mathe´matique, Faculte´s Universitaires Notre Dame de la Paix, Remp
b Dipartimento di Statistica, Universita` Ca’ Foscari, San Giobbe - Cannaregio 873, 3012
c Dipartimento di Scienze Sociali, Cognitive e Quantitative, Universita` di Modena e Reg
d European Center for Living Technology, Calle del Clero 2940, 30124 Venezia, Italy
a r t i c l e i n f o
Article history:
a b s t r a c t
In this paper, we study
journal homepage: wwwll rights reserved.onization in protocell models
tti d
e la Vierge 8, B 5000 Namur, Belgium
nezia, Italy
milia, via Allegri 9, 42100 Reggio Emilia, Italy
eneral protocell models aiming to understand the synchronization
e at ScienceDirect
lsevier.com/locate/yjtbi
etical Biology

Note that synchronization has a further important property,
namely that, even in the case where the kinetic equations for the
whose analyses go beyond the scope of the present work.
ARTICLE IN PRESS
T. Carletti et al. / Journal of Theoretic742GMMs have sub-linear growth terms (Rasmussen et al., 2003;
Munteanu et al., 2006), it leads to exponential growth of the
population of protocells (a straightforward consequence of
constant doubling time)1 and therefore to strictly Darwinian
selection among protocells.
Most of the different protocell architectures which have
been proposed can be divided in two main families, according
to the region of cell space where these key reactions occur.
Some proposals, which have been called surface reaction models
(Serra et al., 2007a)—SRM for short, assume that the key reactions
take place on the surface of the cell membrane, as hypothesized
for the Los Alamos bug (Rasmussen et al., 2003, 2004). Other
architectures exist, for instance the RNA-cell (Oberholzer et al.,
1995; Szostak et al., 2001), where the key reactions develop in the
interior of the vesicle. For this reason we call our model, whose
inspiration has been drawn from this latter case, internal reaction
model—IRM for short.
In this paper, we address the synchronization question for both
proposed architectures, exhibiting a unified analysis; in fact we
are able to prove that working with quantities (of chemicals)
instead of concentrations allows us to map one model on the
other and thus to provide a unified view.
In the case of SRM, we are also able to consider cases
(see Section 3.4) where the ‘‘genetic molecules’’ are actually the
same lipids that compose the ‘‘container’’, allowing us to consider
models close to the GARD—model (Segre´ et al., 1998) or to the one
proposed by Kaneko and Yomo (2002), although in this paper we
limit ourselves to considering only linear interactions. In these
models the (compositional) information is carried by the diversity
of lipids in the vesicle or micelle, thus the synchronization
problem here can be restated in terms of the reproduction of the
whole set of molecules before division occurs, so as to guarantee
the maintenance of information content.
The problem of synchronization has already been studied in
previous works by means of a class of abstract surface reaction
models of protocells (Serra et al., 2007a, b) and it has been
shown that in several cases synchronization is an emergent
property, in the sense that, through successive generations of
protocells, the doubling times of both container and replicators,
tend asymptotically to the same value even if at the beginning
they were different. This was contrasted to earlier models,
like the well–known Chemoton (Carletti and Fanelli, 2007; Ga´nti,
1997; Munteanu and Sole´, 2006), where synchronization was
achieved by ad hoc hypotheses concerning the form of the kinetic
equations.
In models involving a single GMM, synchronization is always
achieved once the growth of the lipid container is linear with
respect to the quantity of the replicator (Serra et al., 2007a). This
result has been generalized to models where the replicator
equation is non-linear or when the growth of the container is
given by a non-linear function of the amount of genetic materia
(Serra et al., 2007b).
In models where more than one GMM coexist in the same
protocell, but limiting the treatment to the case where there is no
direct interaction among them, synchronization was achieved: if
the replicator kinetics is linear, only the fastest replicator
asymptotically survives, while if it is parabolic there is coexistence
of different replicators in the long time limit (Serra et al., 2007a).1 Here we ignore further terms which might limit the growth of the whole
population of protocells, e.g., competition for limited resources or growth in a
limited volume.The treatment of the subject is mostly analytical; we never-
theless present some numerical simulations both to support our
results and to explore cases where the rigorous analysis cannot be
performed.
The paper is organized as follows. In Section 2 we will briefly
introduce the two protocell architectures that somehow inspired
our models: the Los Alamos bug and the RNA–cell. Then we will
introduce our models: the surface reaction model and the internal
reaction model; and finally we will discuss the relevant equations
describing their dynamics. Section 3 will contain a full analysis of
the dynamics of these models and the proof that synchronization
can be achieved, provided some conditions on the involved
coefficients are satisfied. Finally in Section 4 an in-depth
discussion of these conditions and of their physical meaning will
be provided together with some comments on possible further
directions of research.
2. Two protocells models
The aim of this section is to introduce our models describing
two possible architectures for living protocells, inspired by some
current bio-chemical researches. First, for the sake of complete-
ness we will briefly introduce the models and then our approach
will follow.
2.1. Two possible artificial minimal cells
According to Rasmussen et al. (2004), the Los Alamos Bug is a
synthetic organism that integrates three functionalities: a lipid
container, a photo-metabolic system and a hydrophobically
anchored templating polymer that influences metabolic kinetics.
The role of the proto-container is to hold together the other
two key aggregates. Moreover, authors assumed that the interior
lipid phase as well as the water/lipid interface possess very
different physico-chemical properties with respect to bulk water,
in such a way that the high concentrations determined by the
spatial proximity of anchored molecules will enhance the
chemical reactions: i.e., both the lipid phase and the lipid/water
interface act as catalysts.The fact that several different hypotheses lead asymptotically
to synchronization raises the question whether this is a general
property of the SRM or even of a larger class also including the
IRM. In the present paper we therefore explore this wider class of
models taking into account direct interaction, positive and
negative, among the replicators. We consider the case of linear
replication kinetics, finding sufficient conditions to guarantee
synchronization: note however that, since protocell division is
taken into account, the overall model is non-linear, so its analysis
is far from being trivial. We are aware that such assumptions limit
the application of our method to a limited class of models, and
that relaxing them we cannot obtain such analytical results,
nevertheless we support our choice with the following two
reasons. First, we have proved (Serra et al., 2007b) that
synchronization arises under general assumptions of non-linear
coupling between container growth and GMMs or non-linear
kinetics for GMMs replication; second, we stress once again that
we are looking for simplified models able to capture universal
features, neglecting specific, model-dependent details, hence the
linear assumptions are a reasonable starting point. For this same
reason we here neglect higher order phenomena like diffusion and
permeation processes, whose influence can be important but
al Biology 254 (2008) 741–751The assumed GMMs are lipophilic PNA or PNA-like nucleic
acids; this choice has been motivated by the stronger interactions
with the lipid phase of such molecules, thanks to their

hydrophobic backbone; moreover PNA is more plausible than RNA
whe
the
divi
~X ¼
8
ARTICLE IN PRESS
T. Carletti et al. / Journal of Theoreticor DNA in terms of prebiotic synthesis (Nelson et al., 2000). These
polymers, in the single-stranded template form, are located at the
lipid/water surface, exposing the hydrophilic bases to the aqueous
medium while the hydrophobic backbone sinks into the lipid
layer. The system is fed with oligos from the aqueous phase so
that double-stranded templates can be formed—a reaction that is
assumed to be driven by visible light as the primary energy
source—and during this process new lipids are also produced.
Such lipids will be immediately incorporated in the container
triggering its growth, while the double-stranded templates will
eventually dissociate into its two strands, still remaining anchored
to the lipid/water interface. Moreover, the dynamical character-
istics of the reactions are determined by the PNA bases sequences,
hence their role as GMM.
A somehow different approach has been proposed in Oberhol-
zer et al. (1995), Szostak et al. (2001), Mansy et al. (2008),
where the starting point is the design and construction in the
laboratory of a RNA replicase—an RNA molecule that can act both
as a template for the storage and transmission of genetic
information, and as an RNA polymerase that can replicate its
own sequence in the laboratory (Hager et al., 1996; Bartel,
1999; Bartel and Unrau, 1999). This complex must then be
incorporated into some form of compartment, to increase the
replication rate but also because in this way advantageous
mutations can lead to preferential replication: after replications,
mutations and random assortment, some compartments will be
occupied by mutant replicases that can replicate each other
more efficiently, giving them an overall advantage (evolution).
Authors thus hypothesized the existence of a lipid vesicle
membrane, able to self-assemble thanks to the interactions
between the available lipid molecules.
The last step to obtain a living organism is to introduce some
advantage in survival, growth or replication for the membrane
component, directly related to the RNA molecule, for instance a
ribozyme that synthesizes amphipathic lipids pushing the
membrane to grow. Once the container and the genetic material
are coupled the whole protocell will evolve as a whole, and
improved ribozymes—because of favorable mutations—would
have a growth and replication advantage.
2.2. Surface reaction models of protocells
For the sake of completeness in the first part of this section we
recall the model and the main results concerning the surface
reaction models in the case where only one kind of genetic
molecule is present; we refer the interested reader to Serra et al.
(2007a) for more details and demonstrations. Then we generalize
these models to the case of N kinds of GMMs interacting with each
other.
So let us first consider the case where there is a single
replicator, represented by some generic X-molecule in the
protocell lipid phase,2 and let its quantity (i.e., number of moles)
be denoted by X. Let also C be the total quantity of ‘‘container’’
(i.e., lipid membrane forming vesicles or micelles) and V its
volume, which is equal to C=r (where r is the density, assumed to
be constant). Finally, let S denote the surface area, which is a
function of V (S is approximately proportional to V for a
large vesicle with a very thin bilayer membrane, a condition
that will be referred to as the ‘‘thin vesicle case’’, and to V2=3 for a
micelle).
2 This model is invariant with respect to the way in which either amount of C-
molecules or X-molecules are measured; for example, if they were measured as
number of molecules the equations would retain exactly the same form (of course,
the units of the kinetic constants would be different).d~X
dt
¼ Cb1M~X;
dC
dt
¼ Cb1~a  ~X;
>><
>>:
(2.3)
where a ¼ ða1; . . . ;aNÞ is a vector with positive entries denoting
the coupling term between the container growth and each
replicator, while the (constant and real) matrix element Mij
represents the contribution of the XðjÞ-molecule to the growth of
the XðiÞ-molecule. We also assume the matrix M to be invertible, to
avoid redundancy of chemical reactions, and thus to limit the
analysis to the minimal number of independent chemical species.
An important simplification can now be considered: as it wasnon
Eqs.dem
whe
oneolecules. Obviously, all the quantities X must be real and
-negative for i ¼ 1; . . . ;N. The N-dimensional generalization of
(2.1) is thenrepl
~X-mote the total quantity (moles) of N different types of
icating molecules in the protocell lipid phase, called for short
ðiÞdentraightforward. Let
ðXð1Þ;Xð2Þ; . . . ;XðNÞÞ (2.2)T
is svious cycle (perfect halving hypothesis).
he generalization to the case where there are more replicatorsone
preequal daughter protocells and each one will start the next
sion cycle with an initial amount of the X-molecule equal to
-half of the value which it has attained at the end of thegrow
twore Z and a are two positive constants, denoting, respectively,
rate of self-replication of genetic molecules and the container
th. When C reaches a critical threshold, the cell breaks intoWe assume, according to the Labug hypothesis, that the
X-molecule favors the formation of amphiphiles, and that only
the fraction which is near the external surface is effective in doing
so, since precursors are found outside the protocell. We also
assume that the replication of the X-molecule takes place near the
external surface. For example, if the latter is a self-replicating
linear polymer, its replication is accomplished by synthesizing a
complementary chain from free activated monomers. However,
we will develop a treatment which allows greater generality than
pure self-replication.
Following the standard assumptions already used and dis-
cussed in Serra et al. (2007a, b), namely:
(1) spontaneous amphiphile formation is negligible, so that only
the catalyzed term matters;
(2) the precursors (both of amphiphiles and templates) are
buffered;
(3) the surface, S, is proportional to some power of the volume, Vb
(b ranging between 2=3, for a micelle, and 1, for a very thin
vesicle), and therefore, also to the amount of container, Cb;
(4) diffusion is very fast within the lipid phase, so concentrations
can be assumed to be homogeneous;
(5) the protocell breaks into two identical daughter units when its
container size reaches a certain threshold and then halving
the genetic and container molecules between them;
(6) the rate limiting step which may appear in the replicator
kinetic equations does not play a significant role when the
protocell is smaller than the division threshold;
one obtains the following approximate equations which describe
the growth of a protocell between two successive divisions:
dX
dt
¼ ZCb1X and dC
dt
¼ aCb1X, (2.1)
al Biology 254 (2008) 741–751 743onstrated in Serra et al. (2007a), in order to determine
ther there is a synchronization in the asymptotic time limit,
can limit oneself to consider the b ¼ 1 case. The final result

does not depend on b, while of course this parameter affects the
speed with which it is approached; this is essentially a non-linear
rescaling of time, useful to simplify the analysis. With this
simplification, the basic equations (which are valid between two
successive divisions) are then
d~X
dt
¼ M~X;
dC
dt
¼ ~a  X:
8
>><
>>:
(2.4)
We generalize the previous discussion with the following
equations, which refer to the ðkþ 1Þ-th cell division cycle that
starts at time Tk and ends at time Tkþ1:
y
2
¼
Z Tkþ1
Tk
dC
dt
dt and
~Xkþ1 ¼
1
2
~XðTkþ1Þ ¼
1
2
~XðTkþ1; Tk; ~XkÞ, (2.5)
provided Tkþ1 does exist finite, otherwise the division event
will not be defined. Note that in general ~XðTkþ1Þa2~XðTkÞ and
ARTICLE IN PRESS
T. Carletti et al. / Journal of Theoretical Biology 254 (2008) 741–751744As outlined above, we assume that division takes place when
the mass (or equivalently the volume, since density is assumed to
be constant) of the protocell reaches a certain critical size.
Without loss of generality we may then assume that the initial
size of the protocell is one-half of the final value (indeed, if the
size of the very first protocell were different then it would suffice
to consider the evolution from the following generation).
Let us observe here that this assumption, and also the halving
hypothesis of genetic material at division, are not so restrictive as
it could be perceived at first reading. In fact there is nothing
magical about the chosen factor 2, any other factor strictly larger
than one, would give the same results, except of course modifying
the asymptotic values of the amount of ~X-molecules and division
time. As for the second assumption, we are aware that in the
physical case, membrane splitting and divisions of ~X-molecules,
could be better described by introducing some probability
distribution functions of divisions events, with mean value 1=2;
however, this will not change our results except that the
asymptotic values for ~X1 and DT1 will refer to averaged values
of some distribution.
So, starting with an initial quantity of container CðT0Þ ¼ y=2 at
time T0, we assume that once the container size CðtÞ reaches the
critical value y it will divide into two equal protocells of size y=2.
Let DT0 be the time interval needed to double C from this
initial condition, and let T1 ¼ T0 þ DT0 be the time at which the
critical mass y is reached. Clearly DT0 is a function of the initial
quantity of replicators, ~X0. Let us denote ~Xðt; t0; ~X0Þ the quantity of
~X-molecules at time t, i.e., the solution of (2.4) with initial datum
~XðT0Þ ¼ ~X0. When there will be no ambiguity with respect to the
initial datum we are considering, we will use the shorter notation
~XðtÞ. The final value of ~X just before the division is thus
~XðT1; t0; ~X0Þ  ~XðT1Þ. By assumption, each offspring will start the
next cycle with an initial concentration of replicators equal to half
the quantity present at the end of the previous cycle, i.e., in
formula ~X1 ¼ ~XðT1Þ=2, from this point till the next division the
dynamics is determined again by (2.4), let us however observe
that the solution could be different because the initial data have
been set differently. Let us denote the successive doubling time by
T2 ¼ T1 þDT1, thus the third generation will start with an initial
value of genetic material, ~X2 ¼ ~XðT2Þ=2, and so on (see Fig. 1 for a
cartoon describing the construction).Fig. 1. A schematic representation of the construction of the sequence ~Xk . Xðt; tk ;XkÞ de
from an amount Xk at time tk . The division occurs at time tkþ1 and the next protocellDTkþ1aDTk, however we will prove in the next section that these
conditions can be asymptotically approached.
Let us now use the hypothesis that the matrix M is invertible,
so from Eq. (2.4) we get
dC
dt
¼ ~a M1 d
~X
dt
, (2.6)
hence the quantity Q ðtÞ ¼ CðtÞ ~a M1~XðtÞ, is a first integral,
i.e., a quantity constant during each division cycle (the proof is
straightforward, dQ=dt ¼ 0 follows from Eq. (2.6)). Evaluating Q ðtÞ
at the beginning and at the end of the ðkþ 1Þ-th division
we obtain
CðTkÞ ~a M1~XðTkÞ ¼ CðTkþ1Þ ~a M1~XðTkþ1Þ, (2.7)
recalling that C takes the initial value y=2 and the final value y and
using the definition of ~Xk (see Eq. (2.5)) we finally get
y
2
¼ ~a M1ð2~Xkþ1  ~XkÞ. (2.8)
Let us observe that this relation is meaningful only if the container
size increases, otherwise the division event is not even defined
and thus Tkþ1 is not formally defined.
On the other hand, Eq. (2.4) can be explicitly solved during the
ðkþ 1Þ-th division cycle to give
~XðtÞ ¼ ~Xðt; Tk; ~XkÞ ¼ eMðtTkÞ~Xk, (2.9)
hence we get the following relation between the amount of
genetic material between two successive divisions:
~Xkþ1 ¼
~XðTkþ1Þ
2
¼ 1
2
eMDTk~Xk, (2.10)
once again provided Tkþ1 exists.
We postpone the analysis of this model to Section 3, where we
will able to introduce a general framework where the SRM and the
IRM can be put and solved.
Remark 2.1. Let us comment on a simplification which has been
previously used for the SRM model, namely the assumption that
the surface is proportional to a power of the volume. This is
certainly the case for a spherical micelle (with exponent 2=3), but
in the case of a vesicle it holds (with exponent 1) only in the limit
of a very large size and very thin thickness.notes the amount of X-molecule at time t during the continuous growth, starting
cycle will start with an amount of X-molecule given by Xkþ1 ¼ Xðtkþ1; tk;XkÞ=2.

It can be shown that the finite size effects can be taken into
account without modifying our results: synchronization is still
obtained. In fact assuming a generic relation, S ¼ f ðCÞ, between
the surface and the volume, and thus the container size, for
ARTICLE IN PRESS
T. Carletti et al. / Journal of Theoreticsome positive increasing function f, then Eqs. (2.3) have to be
modified into
d~X
dt
¼ f ðCÞ
C
M~X and
dC
dt
¼ f ðCÞ
C
~a  X. (2.11)
But then we can observe that the function given by Eq. (2.6) is still
a first integral and thus the same analysis follows. Another
explanation of this result is that we can ‘‘rescale’’ the time3 by the
positive function C=f ðCÞ and thus identifying Eqs. (2.11) and (2.6).
2.3. Internal reaction models of protocells
In this section we will introduce and analyze a new family of
protocell models inspired by the RNA-cell (Oberholzer et al., 1995;
Szostak et al., 2001). We thus assume the protocell to be a vesicle
made of lipidic material delimiting an inner-water pool, where the
relevant chemical reactions do occur. We also made the simplified
assumption that precursor molecules, of both lipid precursor
and genetic replicators, lying in the environment can enter and
diffuse through the membrane fast enough to consider this step
instantaneous.
We thus have a network of chemical reactions involving the
reproduction of genetic material and the production of membrane
molecules, running in the inner water pool. Let us consider,
however, that the protocell increases its size—the membrane
grows because the produced lipid molecules are supposed to be
instantaneously incorporated into the membrane—hence we have
to account for the volume variations into the ‘‘typical’’ kinetic
equation used. We also assume the following division hypothesis:
once the protocell membrane has doubled its size, with respect to
the initial one, then protocell division occurs, guided by physical
instabilities at the membrane level. At this point two identical
protocells are obtained as offspring under the requirement of
conservation of membrane and genetic molecules.
Let us start with the simplest case, assuming the protocell to be
endowed with only one kind of genetic molecule, whose concen-
tration at time t will be denoted by ½XðtÞ, i.e., the ratio of the
amount, or mass of X-molecules at time t, XðtÞ, divided by
the protocell volume, VðtÞ, enclosed by the membrane. We here
hypothesize the membrane thickness to be very small in such a
way that we can well approximate the internal volume with the
total volume of the protocell, hence ½XðtÞ ¼ XðtÞ=VðtÞ. We also
assume that the membrane surface S is proportional to the quantity
of membrane molecules C and that the volume depends4 on S.
Under the hypotheses of a linear growth for the reproduction
of the genetic material and also that the membrane size
growth is linear in the amount of X-molecules, we get the
following system5:
dC
dt
¼ a½XV ;
d½X
dt
¼ Z½X  ½X
V
dV
dt
;
8
><
>:
(2.12)
3 More precisely let us introduce a new non-linear time t ¼ R t C1ðsÞf ðCðsÞÞds
and let us denote the quantities C and ~X using this new time, respectively, by cðtÞ
and ~xðtÞ, then Eq. (2.11) is formally equivalent to Eq. (2.3).
4 If the protocell were be a sphere then V ¼ ð6
ffiffiffipp Þ1S3=2 / C3=2.
5 Let us observe that assuming the membrane thickness to be constant, one
can express the surface of the protocell as a function of the container C, S ¼ Cd, and
thus the first equation in (2.12) can be rewritten in terms of geometric quantities
as: dS=dt ¼ a0 ½XV .where Z and a are two positive coefficients taking into account,
respectively, the rate of self-replicating of X molecules and the
rate of container growth, which also takes into account the
amount of precursors that we assume to be constant. The last
term in the second equation is due to the varying volume and we
stress once again that we neglect the membrane thickness in
computing the protocell volume.
Let us denote ½Xk the concentration of X at the beginning
of the ðkþ 1Þ-th division and once again DTk ¼ Tkþ1  Tk the
duration of the ðkþ 1Þ-th protocell cycle. The aim is to study
the behavior of ½Xk and DTk as a function of the division
number and show, if any, the existence of some synchronization
mechanism.
Let us start by solving the first Eq. (2.12). Because ½X is never
zero, otherwise everything will stop, the solution describing the
concentration of ½X during the k-th protocell cycle, is given by
½XðtÞ ¼ ½Xk
Vin
VðtÞ
eZðtTkÞ, (2.13)
where Vin ¼ VðTkÞ is the initial volume and ½Xk the concentration
of X-molecules at the beginning of the cycle, while VðtÞ is the
protocell volume at time t. Because we ignore the time evolution
of the volume, and we do not want to impose it, we have to use
the remaining equation of (2.12), that can be easily rewritten
using the second one, as
dC
dt
¼
a
Z
d
dt
ð½XðtÞVðtÞÞ, (2.14)
which implies the existence of a first integral: RðtÞ ¼ CðtÞ
ða=ZÞ½XðtÞVðtÞ. Hence equating RðtÞ at the beginning and the end
of the ðkþ 1Þ-th division we get
CðTkþ1Þ  CðTkÞ ¼
a
Z ð½XðTkþ1ÞVðTkþ1Þ  ½XðTkÞVðTkÞÞ. (2.15)
The hypothesis that volume changes are dictated only by the
surface variations is based on the assumption of turgid vesicle,
which in turn implies a relatively slow membrane growth, in such
a way the volume can adjust to equilibrate the external and
internal pressure. The explicit time variation of the volume is thus
limited to fluctuation around this equilibrium shape that we
neglect assuming it to be small.
We are now interested in determining the concentration of the
X-molecule at the beginning of each division cycle as a function of
the involved quantities. Recalling that, by the halving hypothesis,
each offspring will be endowed with half the number of molecules
produced in the previous cycle, and using moreover the fact that
the concentration is obtained dividing this number by the volume
at the beginning of the cycle, we obtain
½Xk ¼
XðTkÞ
2
1
Vin
¼
½XðTkÞVfin
2
1
Vin
, (2.16)
where Vfin is the volume enclosed by the membrane surface just
before the division. Finally using Eq. (2.15) and calling s the
threshold on the container size at which the protocell division
occurs, we get
s
2
¼ aVinZ ð2½Xkþ1  ½XkÞ.
This relation can be rewritten as
½Xkþ1 ¼
½Xk
2
þ
sZ
4aVin
,
and then iterating back in time we get
al Biology 254 (2008) 741–751 745½Xkþ1 ¼
sZ
4aVin
Xk
j¼0
1
2j
þ
½X0
2kþ1
,

s. First, we prove in S n 3.1 that under the
mptions there always ex a division event, namely
ting at time t ¼ Tk, there ts a positive lapse of tim
 Tk such that CðTk þ ¼ y. This is equivalent
iring the container grow 0. Second in Section 3.2
prove that ~Xk eventually con es to ~w where ~w is prop
1 and ~a  ~w ¼ l1b, and the division intervals c
to a constant value DT ! l1 log a, thus synchronization is
ARTICLE IN PRESS
T. Carletti et al. / Journal of Theoretical Biology 254 (2008) 741–751746which in the limit of infinitely many divisions converges to the
asymptotic value:
½X1 ¼
sZ
2aVin
,
meaning that synchronization is obtained with a constant period
given by
DT1 ¼
1
Z log 2.
Let us now consider the more general case where N kinds of
different genetic molecules are present in the same protocell. Let
us denote by ~½XðtÞ, the vector of their concentrations at time t,
i.e., ~½XðtÞ ¼ ð½X1ðtÞ; . . . ; ½XNðtÞÞ.
Let us again assume linear interactions between the GMMs
and, moreover, that the container growth is also linearly
dependent on the amount of genetic molecules. Thus, we get
the system
dC
dt
¼ ~a  ~½XV ;
d ~½X
dt
¼ M ~½X 
~½X
V
dV
dt
;
8
>><
>>:
(2.17)
where V is the volume enclosed by the surface membrane,6 ~a is a
N-dimensional vector with non-negative entries such that ~a  ~½X ¼Pn
i¼1 ai½Xi represents the contribution to the membrane growth
due to ½X1ðtÞ; . . . ; ½XNðtÞ and the (constant and real) matrix
element Mij represents the contribution of the X
j-molecule to
the growth of the Xi-molecule. Once again we introduced a term
due to the volume variations.
Using the assumption that the matrix M is non-singular, we
can determine the first integral:
RðtÞ ¼ CðtÞ  VðtÞ~a M1 ~½XðtÞ,
once again the following relation between the concentrations at
the beginning of each cycle can be found:
~½Xk ¼
~½XðTkÞVfin
2
1
Vin
,
where Vin and Vfin have been defined previously. Calling s the
surface division threshold, we get
s
2
¼ Vin~a M1ð2 ~½Xkþ1  ~½XkÞ. (2.18)
Once again this relation is meaningful only if the membrane is
increasing, in such a way the existence of Tkþ1 is ensured.
Observe that the second equation (2.17) can be explicitely
solved to give
~½XðtÞ ¼ Vin
VðtÞ
eMðtTkÞ ~½Xk, (2.19)
where the volume variation in time is unknown, and thus
~½Xkþ1 ¼
~½XðTkþ1Þ
2
¼ Vin
2Vfin
eMDTk ~½Xk, (2.20)
provided Tkþ1 do exist.
We thus get a set of equations, cfr. (2.18) and (2.20), governing
the evolution of the IRM, formally equal to the ones for the SRM,
i.e., (2.18) and (2.20). These systems will be analyzed in the
following section.6 Once again we assume that the membrane thickness can be neglected and
assuming it constant we can rewrite the first relation in term of the membrane
surface S.k 1
achieved, in fact with our notations  log a40. Thus the ultimate
fate of both SRM and IRM is to synchronize the replication and
division rates, moreover we can predict, respectively, the
asymptotic amount of genetic material for SRM (and concentra-
tion of GMMs for IRM) in terms of the problem data.
Let us first give a proof in a simplified case which nonetheless
contains all the main ideas and then leave complete proof to
Section 3.3.
3.1. Analysis of the division time
Let us assume M possesses N distinct eigenvalues, l1; . . . ; lN . As
previously stated we hypothesize that l1 is real, positive, withto v7 Because eigenvectors are de
entries we mean that all entries hav
to put them all positive.thatfined up to multiplicative factors, b
e the same sign and thus we can muonverge~
verg ortionalkþ1
requk
s, _C4 , we willT DT Þ to that
star exis e DTk ¼
assu ists for all k
step ectio previous3. Results: synchronization in linear SRM and IRM
The aim of this section is to present our main results, namely to
prove under suitable hypotheses, the emergence of the synchro-
nization phenomenom for both SRM and IRM whose physical
relevance will be discussed in the next section, nevertheless we
anticipate here that all the information is contained in the
chemical constants, i.e., of the matrix M and the vector ~a, and our
method allows us to extract it a priori, avoiding any numerical
simulation. The first step is to rewrite both systems in a compact
unified form. Roughly speaking this can be done because of the
assumption of linear rate equations both for the genetic material
and container growth. Formally let us introduce the following
auxiliary notation:
SRM IRM
~XðtÞ ~XðtÞ ½~XðtÞ
~Xk ~Xk ½~Xk
a 1=2 Vin=ð2VfinÞ
a0 1 Vin=VðtÞ
b y=2 s=ð2VinÞ
In fact (2.8), (2.10), (2.18) and (2.20) can be rewritten as
~Xkþ1 ¼ a eMðtTkÞ~Xk and
b ¼ ~a M1ð2~Xkþ1  ~XkÞ. (3.1)
Let us observe that the Eqs. (2.9) and (2.19) governing
the dynamics during each division cycle can also be cast in a
compact form
~XðtÞ ¼ 2a0 eMðtTkÞ~Xk. (3.2)
The aim of this section is to prove that synchronization holds,
provided M has an eigenvalue l1 real, positive, with algebraic
multiplicity 1 and possessing the largest real part, i.e., l14Rlj for all
the remaining eigenvalues lj. Moreover let ~v1 be the eigenvector
associated to l1 then we assume that ~v1 has positive entries.7
Thanks to the unified approach, we have just introduced, we
can limit ourselves to consider only the SRM case and the IRM
being completely analogous. The proof is thus obtained in twoy positive
ltiply them

ARTICLE IN PRESS
reticalgebraic multiplicity 1 and with the largest real part, i.e., l14Rlj
for all j ¼ 2; . . . ;N. Moreover we assume that its associated
eigenvector ~v1 has positive entries (see Footnote 7).
By standard results of linear algebra the N eigenvectors
of M define a basis of the whole space on which we decompose
the vectors:
~Xk ¼ xð1Þk ~v1 þ    þ x
ðNÞ
k
~vN and
~a ¼ að1Þ~v1 þ    þ aðNÞ~vN , (3.3)
then the solution of (3.2) can be rewritten as
~XðtÞ ¼ 2a0ðel1ðtTkÞxð1Þk ~v1 þ    þ elNðtTkÞx
ðNÞ
k
~vNÞ, (3.4)
and we can easily observe that if t  Tk is large enough the first
term involving the eigenvalue with largest real part will dominate.
Hence the equations for the growth rate of the container can be
rewritten as follows:
_C2a0 el1ðtTkÞxð1Þk að1Þ. (3.5)
Because ~a has positive entries we can assume generically that
að1Þ40, thus if t  Tk is large enough we get the growth of the
container, i.e.,
_C40, (3.6)
which gives the desired result provided also xð1Þk 40 for all k.
We are now able to prove by induction that the division event
always exists. For the first division let us assume that the initial
vector ~X0 has positive projection on ~v1, namely xð1Þ0 40, then by
(3.6) the container will grow and reach the division threshold at
some time T1 ¼ T0 þDT0, during this time ~XðtÞ always has a
positive projection on ~v1, in fact ~XðtÞ ~v1 ¼ 2a0 el1ðtT0Þxð1Þ0 40.
Hence using the halving hypothesis the second cell cycle will start
with a vector of initial conditions ~X1 still with a positive
projection on ~v1, actually: xð1Þ1 ¼ a0 el1ðtT0Þx
ð1Þ
0 40.
We now assume that the protocell undergoes k divisions as
previously described, and we will prove that a further divisionwill
occur. So by hypothesis, we start the kþ 1 cycle with a vector of
initial conditions ~Xkþ1 with a positive projection on ~v1, xð1Þkþ140,
then from (3.6) we get that the protocell container is increasing
and thus at some future time Tkþ1 ¼ Tk þ DTk it will reach
the division threshold, the next cycle will start with a vector
of initial conditions ~Xkþ1 with a positive projection on ~v1,
xð1Þkþ1 ¼ a0 el1DTkx
ð1Þ
k 40. Let us also observe that the division
intervals are strictly positive, DTk40 for all k. Otherwise, an
infinite amount (or concentration) of genetic material will be
required, against any physical reasonability.
This concludes thus the proof of the first part: starting with an
initial amount (or concentration) of GMMs with a positive
projection of the first eigenvector, the division mechanism
never stops.
Let us now prove the claim about the asymptotic value of ~Xk.
3.2. Asymptotic behavior
The aim of this section is to prove that ~Xk eventually converges
to a well defined vector: ~w where ~w 2 spanð~v1Þ and ~a  ~w ¼ l1b,
and that the division intervals converge to a constant value
DTk ! l11 log a, thus providing synchronization.
The second equation of (3.1) can be rewritten as
~a M1ð2~Xkþ1  ~XkÞ ¼ b ¼ ~a 
~w
l1
¼ ~a M1~w, (3.7)
T. Carletti et al. / Journal of Theothus introducing ~dk ¼ ~Xk  ~w, this last relation gives
~a M1ð2~dkþ1 ~dkÞ ¼ 0, (3.8)namely
~a M1~dkþ1 ¼
1
2
~a M1~dk ¼
1
2kþ1
~a M1~d0, (3.9)
so, in the limit of infinitely many divisions we conclude that
~a M1~dkþ1 ! 0. (3.10)
Once again assuming generically that ~a has a positive projection
on ~v1 (and thus also on ~w) and because M is invertible we
conclude that
~dkþ1 ¼ ~Xkþ1  ~w ! 0. (3.11)
To prove that DTk !  log a=l1 let us consider the first equation
of (3.1) and rewrite it using the vector ~dk:
~dkþ1 þ ~w ¼ a eMDTk~dk þ a el1DTk ~w, (3.12)
or, reordering the involved terms:
~dkþ1  aeMDTk~dk ¼ ða el1DTk  1Þ~w, (3.13)
and in the limit of infinitely many divisions, recalling that ~dk ! 0
we obtain
ða el1DTk  1Þ~w ! 0, (3.14)
hence DTk ! 
log a
l1
.
3.3. General case
Let us now rapidly sketch the main changes one has to
consider to deal with the general case where M has only N0oN
eigenvalues, but still assuming that M has an eigenvalue l1 real,
positive, with algebraic multiplicity 1 and with the largest real
part, i.e., l14Rlj for all the remaining eigenvalues lj, and
moreover its corresponding eigenvector has positive entries
(see Footnote 7).
The N0 eigenvectors cannot define a basis of the whole space,
but standard linear algebra results ensure that the set of N0 vectors
can be completed to give a basis using the Jordan vectors. Thus the
decompositions (3.3) must be replaced with:
~Xk ¼ xð1Þk ~v1 þ ~X
?
k and ~a ¼ að1Þ~v1 þ~a?, (3.15)
where ~X?k and ~a? belong to the invariant subspace orthogonal to
~v1. But then using the hypothesis that l1 has the largest real part
and that the Jordan decomposition gives rise to invariant
subspaces it is easy to prove the analogous of (3.6): if t  Tk is
large enough the dynamics is still governed by l1 and thus the
container size increases provided xð1Þk and að1Þ are positive. As the
proof continues similarly, we therefore omit it.
In this section we have thus proved the following result:
Theorem 3.1. Let l1 be the eigenvalue of M with largest real part. If
l1 is positive, with algebraic multiplicity 1 and its associated
eigenvector ~v1 can be chosen with positive components. Then the
SRM model (2.3) exhibits synchronization, provided ~a ~v140 and
~X0 ~v140.
Moreover the asymptotic state is described by
SRM: ~X1 ¼ k~X1k~v1 : k~X1k ¼
l1y
2~a ~v1
and
DT1 ¼
log 2
l1
¼
y log 2
~a  ~X1
.
Let us observe that thanks to our unified approach we also have
the following result:
al Biology 254 (2008) 741–751 747Theorem 3.2. Let l1 be the eigenvalue of M with largest real part.
If l1 is positive, with algebraic multiplicity 1 and its associated

eigenvector~v1 can be chosen with positive components. Then the IRM
(2.12) do exhibit synchronization under the assumptions: ~a ~v140
and ~½X0 ~v1a0.
Moreover the asymptotic state is characterized by:
IRM: ~½X1 ¼ k ~½X1k~v1 : k ~½X1k ¼
l1s
2Vin~a ~v1
and
DT1 ¼
1
l1
log
2Vfin
Vin
¼ s
2Vin~a  ~½X1
log
2Vfin
Vin
,
where Vin, respectively, Vfin denote the protocell volume at the
beginning and at the end of the division cycle.
The physical interpretation of the above results will be provided
in Section 4, here we limit ourselves to emphasize that from the
knowledge of the chemical constants, i.e., the matrix M and the
vector ~a, we can conclude if synchronization will arise or not.
Moreover we are able to define the asymptotic division time and
ARTICLE IN PRESS
T. Carletti et al. / Journal of Theoretic748the amount of each genetic memory molecule. We would also like
to emphasize here that our analysis applies to models involving
chemical reaction where replicators do not interact. More
interesting cases, that can be described by our models, are
represented by replicators that not only positively catalyze each
other’s synthesis but also inhibiting molecules are allowed
provided their ‘‘influence’’ on the chemical network is small
enough.
3.4. Linear GARD—like models
In this last part of the present section, we briefly introduce and
discuss a class of protocells’ models similar to the GARD—model
(Segre´ et al., 1998), where there is no longer a distinction between
genetic material and container, namely the lipids themselves that
form the protocell act also act as information carriers: their
relative amount determine the compositional information
brought by the protocell (see Fig. 2 for a schematic representation
of the model).
We thus assume that lipid precursors are freely available
(without any limitation) in the environment, and that they react
with the lipids at the surface of the protocell to produce new
lipids that are immediately incorporated into the protocell itself,
which thus grows in size. We assume that once the size has
reached some threshold, say twice the initial one, then the
protocell splits into two identical offspring halving the mother
lipid content (i.e., the compositional information). This model can
thus be casted in the SRM class. The amount of the i-th lipid, XðiÞ,
M
ii
M
ij
PrecursorsPrecursors
Fig. 2. A schematic representation of the Linear GARD—like model.evolves in time according to (see Eq. (2.3)):
dXðiÞ
dt
¼
XN
j¼1
MijX
ðjÞCb1, (3.16)
except that now the container size is given by
CðtÞ ¼ Xð1ÞðtÞ þ    þ XðNÞðtÞ ¼ ~a  ~XðtÞ, (3.17)
where we introduced ~a ¼ ð1; . . . ;1Þ, i.e., the total amount of lipid
molecules. We call these models linear GARD—like because the
chemical reactions involved in our model are linear in contrast
with the GARD model where quadratic reactions have also been
considered.
Because the container size is always positive we can ‘‘rescale’’
the time so as to include the factor Cb1 into the new time, hence
the behavior can be described by
d~X
dt
¼ M~X and CðtÞ ¼ ~a  ~XðtÞ. (3.18)
The division and halving hypotheses imply that the amount of
lipids between two successive cycles are related by the follow-
ing law:
y ¼ CðTkþ1Þ ¼ ~a  ~XðTkþ1Þ ¼ 2~a  ~Xkþ1, (3.19)
but
y
2
¼ CðTkÞ ¼ ~a  ~Xk, (3.20)
thus
~a  ~Xkþ1 ¼ ~a  ~Xk. (3.21)
Calling l1 the real, positive eigenvalue of M with largest real
part and ~v1 its associated eigenvector, with positive entries, we
can conclude that following the previous analysis for the SRM
models, the long term dynamics is completely determined by this
eigenvalue and that the asymptotic amount of lipids will have a
positive projection, only on the direction of ~v1, its absolute value
being determined by (3.21). Namely,
~Xk ! ~X1 ¼ k~X1k~v1 such that
k~X1k ¼
y
2~a ~v1
, (3.22)
and
DTk ! DT1 ¼
log 2
l1
. (3.23)
4. Discussion
Let us now analyze the physical conditions ensuring that the
matrix M has a single eigenvalue with largest real part (ELRP for
short) and a corresponding positive eigenvector.
We first discuss the important case where all the matrix
elements are non-negative, i.e., MijX0, for all i; j ¼ 1; . . . ;N. This
implies that there is no negative interference between different
replicators i and j, the only possible alternatives being that either i
favors (e.g., catalyzes) the formation of j or that it does not
influence it in any way. Moreover, we must also require that at
least one of the entries Mij does not vanish, since otherwise there
would be no replication at all. We can therefore apply the
Perron–Frobenius theorem (Lu¨tkepohl, 1996; Minc, 1988), which
states that under the previous assumptions then the eigenvalue
with the largest module is real, positive and unique, and that there
al Biology 254 (2008) 741–751is a non-negative eigenvector belonging to that eigenvalue.
Let us emphasize here that our method enables us also to
predict which replicator molecules will go extinct, in fact once the

our analysis is able to cover several relevant cases, some further
investigations are needed to understand some remaining cases.
Moreover we are aware that non-linear terms may play a key role
and thus further studies are surely necessary to give a compre-
hensive account of the behavior of non-linear systems.
Acknowledgment
Support from the EU FET–PACE project within the 6th Frame-
work Programme under contract FP6–002035 (Programmable
Artificial Cell Evolution) is gratefully acknowledged. Also Elena
Lynch of the European Center for Living Technologies is kindly
acknowledged for the final redaction of the manuscript.
References
Bartel, D.P., 1999. In: Gesteland, R.F., Cech, T.R., Atkins, J.F. (Eds.), The RNA World,
Cold Spring Harbor Laboratory Press, New York, 1999, pp. 143.
Bartel, D.P., Unrau, P.J., 1999. Constructing an RNA world. Trends Cell Biol. 9, 9–13.
Carletti, T., Fanelli, D., 2007. From chemical reactions to evolution: emergence of
species. Europhys. Lett. 77, 18005p1–p6.
Ga´nti, T., 1997. J. Theor. Biol. 187, 583–593.
Golub, G.H., Van Loan, C.F., 1996. Matrix Computations, third ed. Johns Hopkins
University Press, Baltimore, USA.
Hager, H.J., Pollard, J.D., Szostak, J.W., 1996. Ribozymes: aiming at RNA replication
and protein synthesis. Chem. Biol. 3, 717–725.
Kaneko, K., 2006. Life: An Introduction to Complex System Biology. Springer,
The Netherlands.
Kaneko, K., Yomo, T., 2002. On a kinetic origin of heredity: minority control in a
relicating system with mutually catalytic molecules. J. Theor. Biol. 214,
563–576.
Lu¨tkepohl, H., 1996. Handbook of Matrices. Wiley, New York.
Mansy, S.S., et al., 2008. Template–directed synthesis of a genetic polymer in a
model protocell. Nature Lett., doi:10.1038/nature07018.
Minc, H., 1988. Nonnegative Matrices. Wiley, New York.
Munteanu, A., Sole´, R.V., 2006. Phenotypic diversity and chaos in a minimal cell
model. J. Theor. Biol. 240, 434–442.
Munteanu, A. et al., 2006. Generic Darwinian selection in protocell assemblies,
SFI Working Paper 06-09-032.
Nelson, K.E., Levy, M., Miller, S.L., 2000. Peptide nucleic acids rather than RNA
may have been the first genetic molecule. Proc. Natl. Acad. Sci. USA 97 (8),
3868–3871.
Oberholzer, T., et al., 1995. Enzymatic RNA replication in self-reproducing
vesicles: an approach to a minimal cell. Biochem. Biophys. Res. Commun. 207,
250–257.
Rasmussen, S., et al., 2003. Bridging nonliving and living matter. Artif. Life 9,
269–316.
Rasmussen, S., et al., 2004. Proto-organisim kinetics: evolutionary dynamics
of lipid aggregates with genes and metabolism. Origins Life Evol. Biosph. 34,
171–180.
Segre´, D., et al., 1998. Graded autocatalysis replication domain (GARD): kinetic
analysis of self-replication in mutually catalytic sets. Orig. Life Evol. Biosph. 28,
501–514.
Serra, R., Carletti, T., Poli, I., 2007a. Synchronization phenomena in surface-reaction
models of protocells. Artif. Life 13, 123–138.
Serra, R., Carletti, T., Poli, I., 2007b. Surface-reaction models of protocells.
In: Mondaini, R.P., Dila˜o, R. (Eds.), BIOMAT2006. World Scientific, Singapore.
Szostak, D., Bartel, P.B., Luisi, P.L., 2001. Synthesizing life. Nature 409, 387–390.
ARTICLE IN PRESS
T. Carletti et al. / Journal of Theoretical Biology 254 (2008) 741–751 751