Recent Theoretical Approaches to Minimal Artificial Cells

Minimal artificial cells (MACs) are self-assembled chemical systems able to mimic the behavior of living cells at a minimal level, i.e. to exhibit self-maintenance, self-reproduction and the capability of evolution. The bottom-up approach to the construction of MACs is mainly based on the encapsulation of chemical reacting systems inside lipid vesicles, i.e. chemical systems enclosed (compartmentalized) by a double-layered lipid membrane. Several researchers are currently interested in synthesizing such simple cellular models for biotechnological purposes or for investigating origin of life scenarios. Within this context, the properties of lipid vesicles (e.g., their stability, permeability, growth dynamics, potential to host reactions or undergo division processes…) play a central role, in combination with the dynamics of the encapsulated chemical or biochemical networks. Thus, from a theoretical standpoint, it is very important to develop kinetic equations in order to explore first—and specify later—the conditions that allow the robust implementation of these complex chemically reacting systems, as well as their controlled reproduction. Due to being compartmentalized in small volumes, the population of reacting molecules can be very low in terms of the number of molecules and therefore their behavior becomes highly affected by stochastic effects both in the time course of reactions and in occupancy distribution among the vesicle population. In this short review we report our mathematical approaches to model artificial cell systems in this complex scenario by giving a summary of three recent simulations studies on the topic of primitive cell (protocell) systems.


Introduction
The chemical implementation of minimal artificial cells (MACs) from scratch is attracting the interest of a growing number of researchers in the fields of synthetic biology and origin of life studies [1][2][3][4][5][6], who are becoming aware of the potential of microcompartments and lipid vesicle technologies to uncover biologically relevant phenomena, as well as plausible prebiotic processes and evolutionary transitions.A MAC is a self-assembled chemical system that exhibits self-maintenance (i.e.autopoiesis), self-reproduction and the capability of evolution.Self-maintenance represents the capability of an organism to self-sustain itself by transforming chemical compounds available in the environment into building blocks and thus being in a dynamic state of continuous renewal of its constituents.Maturana and Varela formulated in the mid-seventies the theory of autopoiesis in order to define what are the minimal properties of living systems, expressed in terms of their molecular organization [7].Self-maintained chemical structures can also reproduce themselves.This property also stems from their autopoietic organization, but requires that metabolic production of their components overcome their spontaneous degradation.Self-reproducing (proliferating) organisms may exhibit evolution capability in a Darwinian sense as a collective trait due to the transmission, generation by generation, of operational instructions in a coded form.Because reproduction is often imperfect, self-reproducing organisms generate "diversity" and evolve in a Darwinian sense, i.e., owing to the interplay between the generation of genotypic diversity and selection/competition processes in a world of limited resources.While self-reproduction and evolution are typical features of living cellular forms as we know them, it is reasonable to focus on autopoietic self-maintenance for determining the core feature that characterizes minimal "living" entities.
The experimental approach to the construction of MACs is mainly based on lipid vesicles (liposomes) that consist of a closed, spherical, semi-permeable membrane formed by the spontaneous self-assembly of lipid molecules.The membrane is a highly organized molecular bilayer that separates the molecules entrapped inside the vesicle (i.e., in the inner aqueous vesicle core) from the external environment.A variety of chemical and biochemical reactions have been implemented inside liposomes, from RNA synthesis to gene expression, from DNA amplification to lipid synthesis (for a review, see [8]).The latter reaction is particularly important because it allows the growth of vesicles through the enlargement of vesicle membrane.Division may also follow vesicle growth, so that two "daughter" vesicles are obtained from the original parent (i.e., self-reproduction).A schematic classification of artificial chemical structures [9][10][11][12] is tentatively shown in Figure 1, although the nomenclature of these structures are far from being uniformly adopted by the community.In this article, we use the term protocell according to the definition introduced by Ruiz-Mirazo [13], in order to indicate compartmentalized molecular systems (generally lipid vesicles [14], but other compartments have been also used [9,10]) built for understanding the origin of complex biological cells, or for investigating alternative forms of cellular organization [15].As shown in Figure 1, here we will distinguish between truly autopoietic "self-reproducing" systems (capable of growing and dividing by producing the membrane building blocks-the lipids-as well as all catalytic components, for example those present in the aqueous core) and not-autopoietic "self-producing" systems (capable of growing and dividing by producing only membrane components).This distinction is useful because some experimental systems fall in the second category.Relationship between artificial chemical systems aimed at mimicking cellular behavior.Autopoietic systems are all the systems that fulfill the definition given by Maturana and Varela [7] and they can belong to completely different fields, like natural science, sociology or economy.A subset of this collection is represented by autopoietic chemical systems.This subset encompasses reverse [9] and direct [10] micelles (not explicitly shown) along with autopoietic vesicles [11].Depending on the kinetic regime, autopoietic vesicles can self-reproduce, stay in a homeostatic regime or decay [12].MACs can be classified as a subset of autopoietic self-reproducing vesicles since they may also exhibit the capability to evolve.On the other hand, we may call "self-producing" vesicles those vesicles that are able to self-produce at least the membrane, owing to an internal reaction (see the gray set).Whereas autopoietic vesicles are necessarily self-producing, self-producing vesicles are not necessarily autopoietic (their growth being limited to the membrane constituents).According to the Ruiz-Mirazo definition [13] all these chemical systems can be generally indicated as protocells.The three systems discussed in this article are explicitly shown in the diagram (in violet).
In Figure 2 a schematic drawing of a "self-producing" protocell endowed of an internal metabolic cycle that produces the lipid molecules forming the membrane is shown.In this rapidly progressing field, different theoretical approaches have been recently proposed at different levels of complexity, and making the underlying chemistry more or less explicit in the model.Bolton and Wattis focused on early examples of chemical autopoiesis studying the time evolution of the size distribution of autopoietic vesicles [16][17][18] by means of a deterministic kinetic stepwise aggregation model: The Becker-Döring model.In this approach a vesicle is described only in terms of its aggregation number or its size, i.e. the number of lipids that form the membrane, and by using a coarse grain reduction technique, these authors were able to reduce the entire vesicles population to a few aggregates groups of various sizes each described using a single time variable.

Figure 2:
Schematic drawing of a protocell capable to synthesize membrane constituents ("self-producing" system).The internal metabolic cycle produces fresh lipids L by converting lipid precursors X present in the external environment.The vesicle can grow and divide only if the lipid production is faster than the lipid degradation: growth regime, otherwise, if these processes are balanced, a homeostatic regime takes place or, if lipid degradation is the most efficient process, vesicles decay in time.

Membrane
Aqueous Core Transport Process

Internal Metabolism
On the other hand, biophysical studies on vesicle self-reproduction have been reported by Bozic and Svetina [19,20] and have focused on morphological shape transitions of a single self-reproducing vesicle.These authors analyzed the physical conditions under which reproduction would take place under the assumption that chemical reactions provoke an exponential increase in the total membrane area of the vesicle.By minimizing the bending energy of the membrane, they were able to derive a deterministic equation for division in terms of the different physical magnitudes involved (permeability, spontaneous curvature, bending constant).With a similar approach, Fanelli and McKane also developed alternative models [21].In contrast, Surovtsev et al. [22,23] proposed a much more explicit chemical reaction network and studied how internal chemical oscillations are critical to induce adequate shape transformations and membrane division.Making use again of the principle of minimal bending energy, they showed, in this case, how system-controlled fission is nearly impossible to accomplish without a rather sophisticated mechanism (e.g., A contractile ring to ensure the effective narrowing of the division neck).More elaborate models have been proposed by Macia and Solé [24] who explored how chemical/metabolic symmetry-breaking concentration patterns can generate osmotic differences across the global volume of the vesicle, inducing changes in its morphology and, eventually, a division process.The types of reaction networks they studied were reasonably simple and embedded in a spatially explicit 2D model of membrane growth and deformation [25].
We have been involved in the field of synthetic cell modelling with mathematical models implemented by a specialized Monte Carlo platform called ENVIRONMENT [26] which tries to couple the membrane energy state with the internal metabolism of a single protocell or of population of protocells.The aim of this paper is to summarize our recent studies on MACs, and it is not intended to be a comprehensive review on the entire topic.In particular, we aim at introducing the mathematical framework we used to describe the time behaviour of reacting protocells in terms of the deterministic versus the stochastic approach [27].We also would like to focus on a recently emerging theme, namely the interplay between internalized reactions, vesicle growth and self-reproduction.
In our approach, protocells are seen as compartmentalized chemically reacting systems that can change their size and morphology during their lifetime.These variations, which are well evident in terms of reaction volume and membrane area, can highly influence the protocell time behaviour, affecting for instance the integral reaction rates and transport processes.Moreover, since the rate of change of the aqueous core volume and the lipid membrane surface are not directly coupled, a rapid increase of the core volume can bring a protocell into an unstable state and eventually to an osmotic crisis, while conversely a rapid growth of the membrane area can bring vesicles to a deflated state which will eventually result in protocell division.In order to grasp this complex phenomenology, we have introduced a simplified model of reacting vesicles that can be reduced to deterministic equations suitable for the description of the protocells' time-averaged behaviour.Furthermore, the size range of a population of lipid vesicle protocells can be very broad (from tens of nanometers to tens of micrometers).The distribution of biomolecules encapsulated inside vesicles can be largely affected by randomness, giving rise to a population of individual compartments very different from one another in their chemical composition.In order to take into account the role of stochasticity in the protocells' time behaviour, a novel theoretical approach is necessary.This will allow the study of an entire population of protocells instead of a single compartment, like in the deterministic approach.In turn, such an approach might be useful to test and then overcome the main assumption of the deterministic approach, namely that the time behaviours of all the protocells in the population are, on average, statistically equivalent.
This theoretical review is structured as follows: in Section 2 the deterministic and stochastic approaches will be presented and our in silico reacting vesicle model will be used to elucidate the main difference between them, while in Section 3 some examples of protocell modelling will be presented, showing how random fluctuations can bring the system very far from the deterministic forecast.Finally some conclusions will be drawn and an outlook of future work presented.

Mathematical Approaches
In this section, we will introduce and discuss a lipid vesicle model that allows the description of reacting vesicles.The corresponding mathematical equations will be presented according to deterministic and stochastic approaches.In both approaches we will adopt the "continuous stirred tank reactor" assumption for all microcompartments, which is commonly adopted in biological modelling, from environmental analysises to pharmacokinetics studies.

In Silico Reacting Vesicles
According to the schematic drawing shown in Figure 2 self-producing vesicles are described as compartmentalized systems made up of two different homogeneous domains: the membrane and the aqueous core [28].Lipids can be exchanged between the membrane and aqueous core and between the membrane and the external aqueous environment.Transport processes can also occur, allowing for the exchange of molecules from the external environment to the internal water pool.The vesicle membrane surface S µ can be determined as it follows: where the contributions of different lipids present in the bilayered membrane is calculated by the product of hydrophilic head area  i, times the number of lipids i n  and neglecting the membrane thickness ~4 nm (this assumption works better for vesicles with diameter larger than 30 nm).In the rest of this paper we will deal with vesicles whose membranes are made of a single lipid so that the previous equation simplifies to In the presence of osmotic pressure unbalance, the internal aqueous volume V C is considered to be affected by a water flux, so that the variations of the core volume and the membrane surface are not necessarily directly coupled.The vesicle state can be described by the reduced surface : defined as the ratio between the actual membrane surface S  and the area of a sphere with the actual vesicle core volume V c , (Figure 3).Therefore,  equals 1.0 for spherical vesicles, while it is less than 1.0 or greater than 1.0 for inflated or deflated vesicles, respectively.In fact, a flux of water can take place across the lipid membrane driven by an osmotic pressure unbalance.An osmotic rupture takes place when the internal volume increases until eventually the membrane rupture when  < (1−), with  being the osmotic tolerance.
Parameter  is experimentally derived, and determined by measuring the maximum increase ∆C max in the internal concentration C that vesicles can support without breaking [29]: On the other hand, in the present model, deflated vesicles are assumed to divide when the membrane surface is large enough to form two twin spherical daughters: 3

 
. This event has been observed under some experimental conditions [30].However, the dynamics of a deflated vesicle are not at all trivial, and modelling the division process will become more accurate as more experimental data become available.(a)  > 3/2 osmotically-stressed growth, i.e. the volume increases faster, then it will reach an elastic tension condition and, above the limit of elasticity of the membrane, this will lead the vesicle to osmotic burst ( < 1−); (b)  = 3/2 continuous spherical growth, i.e. a spherical vesicle will increase its size without any change of shape ( = 1); (c)  < 3/2 reproductive growth, i.e. the surface increases faster than the two previous cases, the growing vesicle will become deflated, assuming a non-spherical shape (ellipsoidal, elongated or, generally speaking, a prolate shape) and the energy of the membrane will be higher due to a bending tension.
Since the aqueous core volume V C and the membrane surface S µ may follow independent time trends, in order to describe the various possible behaviours of the system, it is convenient to introduce the growth control coefficient  [31]: This dimensionless observable is defined as the ratio between the relative rate of volume variation, and the relative rate of surface variation.In presence of an endogenous (biosynthetic or protometabolic) production of lipid dS > 0 thanks to the spontaneous uptake of fresh lipid by the membrane and, in these conditions, >0 indicates a real growth regime.Therefore, just by applying some straightforward geometry rules for a generic growing sphere of radius R: d(lnV)/3 = d(lnS)/2 = d(lnR), three different scenarios among all possible growth regimes may be distinguished (Figure 4).

Deterministic Modelling
If in the vesicle aqueous core, N species X i (i = 1,2…N) react according to  chemical elementary reactions: then the average time evolution of the reacting vesicles can be described in terms of the average number of molecules i x of internal aqueous species X i (i=1,2…N, i≠L), the lipid molecules L x in the aqueous core, the lipid molecules in membrane L x  and core volume V C by solving the following ordinary differential equation set (ODES) [27]: where N A is Avogadro's number, k in and k out are the kinetic constants of the lipid uptake and release to/from the membrane and the external and internal aqueous solutions, while i  is the membrane permeability and [X i ] Ex and x i /(V C N A ) are the external and internal concentrations of i-species, respectively.The term v  is the reaction rate of the -th reaction ( = 1,2…) of the internal metabolism given by the mass action law: where k  is the kinetic constant.It is worthwhile to stress that the ODE set in Equation ( 6) has been written for the case of aggregates formed by a single type of lipid X L .Others simplifying assumptions include: (i) the molecular diffusion of species both within the internal core volume and in the external environment are neglected; and (ii) the external concentrations [ ]   i Ex X are considered constant in time (i.e., the environment is considered as an infinite source of external compounds and a sink to waste).
Going into details, the molecular number of each aqueous species x i changes in time due to the internal metabolic reaction and to the transport process from the outside.The transport across the membrane is driven by a concentration gradient as shown by the following scheme: where is the internal concentration of the ith species .Instead, the rate of change of the lipid number in the core, L dx dt , takes into account the exchange between the aqueous internal phase and the membrane described as follows: where is the internal concentration of the lipid.Moreover, the lipid exchange towards the outside is not explicitly considered in the rate equation of the membrane lipids L dx dt  since the external lipid concentration is assumed to be constantly equal to the equilibrium value . The last equation in the ODE set Equation describes the core volume change rate due to a flux of water driven by the difference of the total osmolite concentration (i.e. an osmotic pressure unbalance, being  aq the water permeability and  aq the water molecular volume).It is important to remark that the deterministic approach gives the time evolution of the vesicles solution as the average time course calculated over the vesicles population, so that although i x (I = 1,2…N) and L x  are not positive integer numbers but positive real values, nevertheless they represent amounts of molecules.Therefore the vesicle state is represented by the array X = (x 1 , x 2 , …, x N ) T and the core volume V C .When the condition for division is satisfied ( ), then the vesicle divides in two twin daughters with volume equal to V C /2 and all the elements of the state array are accordingly divided by 2.

Stochastic Simulations
The stochastic kinetic approach explicitly takes into account the discrete nature of molecule numbers and the intrinsic randomness of reacting events.Therefore, the state of a reacting vesicle is defined by an array of integer molecular numbers n i : N = (n 1 , n 2 ,… n N , L n  ) T and the core volume V C .Moreover, for each elementary reacting event a propensity density probability a  (N) is introduced instead of the deterministic reaction rate so that a  (N)dt gives the probability that the -th reaction will take place in the next infinitesimal time interval dt [27]: while the propensity density probabilities for transport processes and lipid exchange can be predicted according to Equations ( 5) and ( 6) [26].The stochastic time evolution of a well stirred chemically reacting system can be then obtained by solving the following master equation (ME) [27]: that expresses the change rate of the Markov density function P(N, t| N 0 , t 0 ), i.e. the density probability to find the system in the state N in the time interval [t, t + dt) given the system in the state N 0 at time t 0 .∆N  is the jump array, that is the stoichiometric variation of the number of molecules due to the -th reaction.By solving analytically the ME, the average time behaviour of the reacting system can be obtained, along with displacements from the average species time course due to random fluctuations.These fluctuations can bring the system towards regimes unpredictable by the deterministic approach [27].ME is very difficult to solve analytically, but it can be simulated exactly by the well know Monte Carlo Direct Method introduced by Gillespie [32].Based on this method we developed a software platform [26] suitable for simulating the stochastic time evolution of a collection of reacting vesicles, assuming that spatial diffusion processes can be neglected and the concentration gradients take place only across the lipid membranes.Moreover, this program allows one to study the case of vesicle self-reproduction, since it is able to follow a collection of reacting compartments that increases in number (see [26], and [33][34][35] for further details).We wish to remark here that this software is also suitable to study the influence of extrinsic stocasticity.In fact, reacting molecules can be distributed randomly among compartments at the starting time or between daughters at the division time ( ), simulating how this source of randomness affects the system time behaviour.

Results and Discussion
In this section three examples of protocells will be reviewed and discussed.Firstly, we will present a model of autopoietic oleic acid vesicles in a homeostatic regime.This study, based on experimental data [12], focuses on a minimal model that can describe the observed behaviour by also taking into account the role of random fluctuations.Then, we will derive a general formula suitable for predicting when the stationary self-reproduction of protocells can take place according to our assumptions.This formula will be applied to the simple case of enzymatic self-producing vesicle, a short name to indicate those lipid-synthesizing vesicles that are capable of growth and division due to the enzyme-catalysed conversion of a molecular precursor in the membrane-forming molecule L. Note that in these vesicles the enzyme itself is not produced by any internal mechanism.Despite its apparent simplicity, no experimental data are currently available for such systems.Finally, a more complex case will be presented: RNA-based protocells ("ribocells") which are self-reproducing protocells based on ribozyme catalysts.By means of numerical simulation, it will be checked whether such a hypothetical minimal cellular model [36] is robust enough to be implemented in a test tube.

Autopoietic Vesicles in Homeostatic Regime
Autopoiesis, as developed by Maturana and Varela in the seventies [7], is a theoretical description of the 'blue print' of cellular life.It poses as a main feature the self-maintenance of the cell, as due to a process of components' self-generation from within the cellular boundary-a boundary which is itself one of the products.From the chemical point of view, the fertility of autopoiesis theory allowed the design and the experimental achievement of some autopoietic chemical systems all based on surfactant self-assembling structures, such as micelles [10], reverse micelles [9] and vesicles [11].In Figure 5, the schematic representation of an autopoietic vesicle is shown, along with the kinetic conditions for different regimes in the time course of total surfactant concentration.The occurrence of these regimes depends on the rates of amphiphile production v G and decay v D , respectively.Zepik et al. [12] have reported experimental work that confirms the expected behaviour.In particular, the chemical system consists of a solution of oleic acid/oleate vesicles (S i , i being the aggregation number), buffered at pH 8.8, fed with a surfactant precursor and with a reactant capable of destroying oleic acid.The surfactant precursor (P) is oleic anhydride, a hydrophobic substrate rapidly taken up by oleate vesicles at their membranous interface.Owing to the high pH value, P is converted to oleate by alkaline hydrolysis that takes place on the membrane of vesicles.Oleate vesicles also undergo a decay process due to the simultaneous transformation of oleate molecules into 9,10-dihydroxystearate (W) by osmium tetroxide/potassium ferrocyanide oxidation (Y).The dihydroxylated compound W does not form vesicles; therefore, the consequence of the latter conversion is a stepwise vesicle collapse (death).Due to the two competitive reactions, the overall oleate concentration increases, remains approximately constant, or decreases, depending on the magnitude of the P and Y flux rates [12].In order to reproduce the experimental observed behaviour we proposed the simple mechanism reported on the right side of Figure 5 and we were able to obtain the time course of the overall oleic acid concentration [35]: By performing simulations of autopoietic vesicle populations (each characterized by a specific vesicle size), it has been demonstrated that stochasticity selects vesicles with aggregation numbers in the range 10 3 -10 4 (Figure 6a).This effect can be ascribed to the presence of random fluctuations in the growth and decay specific rates.In real (chemical) reacting systems these fluctuations are due to the intrinsic stochasticity of reacting events, but they can also be enlarged by natural changes of physical parameters such as temperature, molecular fluxes, etc.In fact, stochastic simulations starting from a single aggregate have shown how random fluctuations at the steady state can drive the evolution of the aggregate towards growing or a decreasing in size (Figure 6b).Therefore, when autopoietic vesicles of different sizes are present in a system in stationary conditions, fluctuations can act as a selection rule that leads to the perpetuation of those aggregates large enough to overcome large deviations.In conclusion, stochastic simulations have shown that, in this landscape, random and driven fluctuations can represent the driving force for aggregates evolution, growth or decay, and at the same time they can act as a selection rule for the fittest, i.e. the most robust, aggregates in a prebiotic environment.Reproduced with kind permission of the authors [35] and under the conditions of IOP Publishing.

Protocell Stationary Self-Reproduction
In a recent work, we derived a phenomenological law that predicts when a stationary self-reproduction takes place for minimal autopoietic vesicles.By "stationary self-reproduction" we mean a dynamic regime where the conditions for division are reached repeatedly after a constant, characteristic period of time, giving as a result two vesicles or protocells with the same (initial) size, lifetime and metabolite concentration profile as the progenitor.
In terms of the growth control coefficient it is easy to prove that the steady condition takes place when  = 1.Then, two general expressions for the temporal behaviour of the protocell surface and the protocell core volume have been independently derived [31] and an explicit relationship among different molecular and kinetic parameters (e.g., reaction rates v  , permeability coefficients  i , metabolite concentrations [X] i ) have been analytically derived for the protocell stationary reproduction: where v L is the rate of lipid production, C C is the total internal concentration and ∆m  is total variation of the number of molecules due to the -th reaction: , , ( ) Equation ( 13) shows the deterministic condition for a stationary reproduction regime that results from the osmotic synchronization between membrane and core volume growth, i.e.: A spontaneous "self-regulation" driven by the osmotic balance across the protocell lipid bilayer.Equation ( 13) links metabolic kinetic constants and membrane permeabilities with the external and internal concentrations of the system constituents.Therefore, it represents a constraint for the possible sizes and division periods of stationary self-reproducing protocells.We have applied the general Equation ( 13) to the simplest case of a self-producing enzymatic vesicle (SPEV) represented in Figure 7. SPEV is a hypothetical protocell model where the production of lipid S takes place through the chemical transformation of a precursor molecule P, assumed to occur only in the presence of an additional compound E (i.e. an enzyme) encapsulated in the core volume.The S production generates also the waste W so that ∆m = 1 and the osmotic synchronization can in principle takes place.Moreover, W is accumulated in the core volume since it is assumed not to be transported across the membrane, i.e.  W = 0.It is worthwhile to note that this model is very close to some experimental approaches based on giant vesicles that produce internally (with the help of a synthetic catalyst) the main membrane component and eventually undergo growth and division [37].SPEV is not a full-fledged autopoietic vesicle since the catalytic specie E is not synthesized by the internal metabolism.Therefore after each vesicle division the number of E molecules will decrease until just one copy of this molecule is present in the internal core.As a consequence, whenever a division occurs only one of the two daughter vesicles will be able to encapsulate the catalyst molecule and, therefore, will keep the potential to continue growing, producing S and, eventually, dividing in two.The vesicle that contains that single molecule E, by default, will be taken as the mother vesicle, whereas the daughter (and all possible granddaughters) will be "sterile" vesicles.Thus, by using Equation ( 13) it was possible to predict [31] for the mother SPEV, i.e. the vesicle that keeps the only E molecule, its stationary radius R  : and the division time ∆t  : where C C is the overall internal osmolite concentration, [P] Ex and  X are the external concentration and the membrane permeability of the lipid precursor respectively, while k is the kinetic constant of the lipid production: Figure 8 shows in the left plot the core volume time trend for the first seven generations, i.e. vesicle divisions, obtained both by ODE set integration (black line) and by stochastic simulations (gray lines).Vertical gray dotted lines represent the time of division that takes place when the reduced surface satisfied the splitting conditions: .Generation by generation the mother protocell tends to the stationary growth and division as illustrated by the left plot where the values of the vesicle volume before 2V  and after the division  have been calculated with Equation ( 15).In the right plot it is reported the division time ∆t g against the generation number, showing that generation by generation it tends to ∆t  as predicted theoretically.An important aspect to remark is that Equation ( 13), strictly speaking, only captures the condition for stationary reproduction in the sense of a global synchronization process between membrane and volume growth.In other words, it does not guarantee that when a vesicle reaches the division threshold the number of each internal constituent gets effectively doubled (with regard to its initial state in the protocell cycle).This becomes manifest in the case of SPEV, where the single enzyme/catalyst present in the mother is not doubled and, therefore cannot be transferred but to one of the offspring vesicles (i.e.: The only one that will remain fertile).Therefore, Equation ( 13) states a necessary but not sufficient condition for reliable reproduction of proto-cellular systems.In a more complex scenario, which will be introduced in the next section, the metabolic reaction network includes also the synthesis of the enzymatic/catalytic compound and the full reproduction of the entire protocell can be achieved.But the synchronization among lipid production, enzyme duplication and membrane division only emerges if the metabolic pathway(s) led to effective internal chemical synthesis (i.e.∆m>0), since the mechanism that drives the synchronization is the osmotic balance across the lipid bilayer.Moreover in this complex scenario stochastic fluctuations also affect the time behaviour of each single protocell significantly more, since they can strongly influence the distribution of enzymatic species between the two daughter vesicles.

A Minimal Cell Model: The Ribocell
The so-called ribocell (RNA-based cell) is a theoretical minimal cell model based on a self-replicating minimum RNA genome coupled with a self-reproducing lipid vesicle compartment that has been recently hypothesized [36].This model supposes the existence of two ribozymes, one (the lipid synthase R L ) able to catalyse the conversion of molecular precursors (P) into lipids (S) and the second (the polymerase R P ) able to replicate RNA strands by a template driven elongation.Therefore, in an environment rich in both lipid precursors (P) and activated nucleotides (NTPs), the ribocell can grow and self-reproduce if both processes, i.e. genome self-replication and membrane production are somehow synchronized, leading to a sustainable and recursive growth and division dynamics.Recently we have explored the feasibility of this hypothetical minimal cell [33] by determining the best external conditions to observe synchronization between genome self-replication and vesicle membrane reproduction, making use of a deterministic kinetic analysis, while the ribocell robustness to random fluctuations has been tested by stochastic simulations.The proposed metabolic mechanism is reported in Figure 9.Both pairs of RNA strands reversibly associate (A).The association is shifted towards the dimer formation and is strongly dependent on temperature.The replication of every RNA strand is catalysed by the polymerase R P according to the steps in bracket (B).This process is described as a catalytic template-directed addition of mononucleotides with high fidelity and processivity.It starts with R P binding any of the monomeric templates T (T=R P , c R P , R L and c R L ) to form the complex R@T.This complex will then initiate the polymerization of the conjugate strand c T, by iteratively binding NTPs by means of specific Watson-Crick interactions, joining the nucleotides together to form the novel RNA strand c T, and releasing the by-product W. When the strand c T has been completely formed, the polymerase ribozyme releases the new dimer.Finally, the ribozyme R L catalyses the conversion of the precursor P into the lipid S (C).All the kinetic constants have been estimated by experimental values reported in literature and are listed in Table 1 along with references [38][39][40][41][42][43].
Thanks to a deterministic analysis [44,45], we showed that if the kinetic constant for lipid formation k L is in the range: 1.710 3 s −1 M −1 ≤ k L ≤ 1.710 5 s −1 M −1 then synchronization between vesicle reproduction and genome replication can spontaneously emerge under the model assumptions and kinetic parameters reported in Table 1.Deterministic calculations were performed for two ribozymes 20 bases long and showed that the ribocell reaches a stationary growth and division regime ( = 1), where the cell size remains constant after each division along with the amount of genetic material.Although the observed cell life time stabilizes after the first 10 generations, it remains very high (about 80 days for all the k L values in the synchronization range), making the ribocell very hard to implement and study experimentally.Therefore, we investigated the robustness of the stationary growth and division regime of the ribocell in terms of the external substrate concentrations, vesicle size and initial ribozyme amount in order to define optimal external conditions for ribocell self-reproduction [33].The influence of ribozyme length was also explored in the optimal external conditions by varying strand size from 20 to 200 bases in length and keeping all the other kinetic parameters constant.An RNA length of 20 bases is in fact the minimum length required to observe a folded RNA structures, i.e. a structure that can reasonably exhibit catalytic action.On the other hand, entities of about 200 nucleotides have been suggested as plausible ancient proto-ribosomes [46] even though, more recently, smaller subunits of 60 nucleotides have also been considered as plausible candidates [47].Therefore, setting the RNA length of both ribozyme pairs to 100 bases, the deterministic analysis shows that starting from external concentrations [NTP] Ex = [P] Ex = 10 −2 M at the stationary regime the ribocell radius is 113.0 nm and the division time reduces to 68.2 days.The total number of RNA strands is 258 and the genome composition is quite uniform 25.2% (R L ), 25.2% ( c R L ) 25.6% (R P ), 24.0% ( c R P ).The stationary division regime can be reached starting from initial genome composition ranging from 1 to 100 dimers of R c R L and R c R P .In Figure 10, the deterministic time behaviour of the ribocell in optimal external conditions is reported.

Figure 9.
The ribocell model: a schematic drawing on the left, and the internal metabolic mechanism in details on the right.R L and R P are the lipid synthase and the polymerase ribozymes respectively, while c R L and c R P are their complementary filaments and R c R L and R c R L are the double strands formed according to the reversible dimerization (A).The scheme (B) is the mechanism for the template driven synthesis of RNA strands where R P @T is the complex polymerase-template and T is the template thah can be R L , c R L , R T or c R T .The reaction (C) is the conversion of the lipid precursor P into lipid S by producing a waste molecule W.
Finally, the dependence of t 25 (division time after 25 generations) on the kinetic constants for RNA dimer formation k TT and dissociation k T has been also studied.It turns out that the ribocell life cycle at stationary regimes does not depend explicitly on the absolute values of kinetic constants k SS and k S but on their ratio: k SS /k S , i.e., on the thermodynamic constant of RNA dimerization.The more thermodynamically stable the RNA dimers, the longer it takes to observe ribocell self-reproduction.
For instance, if k SS /k S is decreased by two orders of magnitude, the ribocell life time reduces from 68.2 days to 6.4-11.8days.
* k L is 10 5 times larger than the value of the splicing reaction, catalyzed by the hammerhead ribozyme.Stochastic simulations have also been performed in order to test the robustness of ribocells based on 100-base length ribozymes in optimal external conditions, with the aim of elucidating the role of intrinsic and extrinsic stochasticity on the time behaviour of a ribocells population.Simulations were executed by means of the parallel version of ENVIRONMENT [26], running 32 statistically equivalent simulations of a 10-ribocell solution on different CPUs.Therefore, the outcomes were obtained as averages from a population of 320 vesicles.Kinetic parameters used for simulations are those reported in Table 1.At each cell division, only one of the two offspring was kept while the other was discarded in order to reduce computation time, thus keeping the number of monitored vesicles constant.This is in agreement with the assumption that the external concentrations of all substrates are fixed due to an incoming flux of material, i.e. the substrates can never be exhausted.The simulation outcomes are reported on the left of Figure 11 where the composition of the ribocell population is reported against time.In fact, during simulations at each division the genetic material is randomly distributed between the daughters.If the amount of genetic material is very low, then this can result in a separation of R P from the other RNA strands.In fact, the ribocell must contain a minimum genetic kit of three RNA filaments in order to be capable of self-replicating its entire genome: one R P that catalyzes the RNA base pair transcription, one (R L or c R L ) and one (R P or c R P ) that work as templates for the transcription.Moreover, since R L is necessary to catalyze lipid precursor conversion, the optimal minimal 3-ribozyme "kit" must be made up of 2R P and one R L .This minimum kit should be at least doubled before cell division, in order to have a chance that both daughters continue to be active.Therefore, if a random distribution of RNA filaments takes place after vesicle division, ribozyme segregation between the two daughters might occur.Different scenarios can be envisaged: death by segregation is reached if vesicles are produced without any ribozymes (empty vesicles), or are produced containing one lone R P , or many filaments of c R P and/or c R L (inert vesicles).Vesicles that encapsulate R L strands are self-producing as defined in his work: they are able to synthesize lipids and then can grow and divide producing in turn self-producing and/or empty vesicles.On the other hand, vesicles containing more than one molecule of R P or both R P and c R P filaments are able to self-replicate this reduced genome (self-replicating genome vesicles) but they cannot self-reproduce the membrane.So their fate is an osmotic burst due to an unbalanced increase in waste concentration.Finally, a reduced version of the Ribocell consists in a lipid aggregate that contains one R P filament and R L / c R L strands.As a consequence of this, reduced ribocells are able to replicate the R L / c R L genetic material, and at the same time are able to synthesize lipids.Therefore, they can grow and divide, producing in turn at least one reduced ribocell and/or self-producing, inert and empty vesicle.
On the right of Figure 11 a schematic drawing of the different types of protocells is reported.At the end of the simulation, the composition of the protocell population consist of low percentages of still fully functioning ribocells (6.7%) while the most populated fractions are those of empty (40.0%) self-producing (33.3%) and broken (20.0%) vesicles, respectively.Reduced ribocells are only present in the first generations since they very soon decay into self-producing and empty vesicles.Inert vesicles, i.e. vesicles entrapping free chains of c R P and/or c R L or a single R P , are not formed and this can be ascribed to the high stability of RNA dimers and complexes (i.e. the chance of finding free RNA monomers at the time of vesicle division is extremely improbable).Indeed, the stochastic time trend presents a very irregular time behaviour compared to the deterministic one that describes a highly synchronized oscillating regime of growth and division.In contrast, stochastic simulations highlight the alternation of dormant phases, where the reduced surface remains practically constant (both the core volume and the membrane surface being constant, data not shown), to very active steps where ribocell growth takes place very fast, leading to a division event.The fast growth and division step corresponds to the presence of a free R L chain in the vesicle core while, in the dormant phase, ribozymes are all coupled in the form of dimers or complexes.As a consequence, self-producing vesicles with a genome made up only of R L monomers can reproduce very efficiently since no dormant phase can occur, given that the formation of R c R L dimers is impossible.These protocells can then self-produce very efficiently, with a ∆t less than one day.[33] with kind permission of the author.Different protocells can emerges due to ribozymes segregation among the daughter protocells at the moment of mother protocell division: vesicles without any ribozymes (empty vesicles) or containing one lone R P or many filaments of c R P and/or c R L (inert vesicles), vesicles encapsulating R L (self-producing), vesicles containing more than one molecule of R P or both R P and c R P filaments (self-replicating genome), aggregates containing one R P filament and R L / c R L strands (reduced ribocells).
In conclusion, the simulation outcomes show that ribocells are not robust enough to survive to random fluctuations.In fact only about the 5%-7% of the initial population survive as full ribocells after 15-25 generations and on a longer time frame they too are destined for extinction.Furthermore, the time course of each single protocell is also greatly influenced by intrinsic stochasticity in particular by the time fluctuations of the RNA dimer dissociation.In fact, when all the RNA strands are associated in dimers, protocells remain in a dormant phase, whereas free R L monomers induce fast growth and division steps and free R P filaments cause the fast RNA replication without changing the vesicle size appreciably.Therefore these two processes are synchronized only by chance and this also represents a reason of weakness of this model protocell.To improve this behaviour, one possibility could be to increase the working temperature in order to facilitate the dimers dissociation and increase the concentration of the free catalysers R L and R P .For further details the reader is addressed to [33,[44][45].

Conclusions
In this short review we have presented and discussed some aspects of theoretical modelling in micro-compartmentalized systems, in particular by discussing our recent research on self-reproducing protocells.The occurrence of compartmentalized synthetic reactions coupled with the membrane dynamics (in terms of growth and division) plays a major role in determining the evolution of these kinds of chemically reacting systems.
We started by comparing the deterministic and stochastic approaches for modelling such systems, and reviewed how these methodologies can be applied to describe: (1) homeostatic autopoietic systems, (2) the stationary conditions for protocell self-reproduction, and (3) the more complex case of the "ribocell".In all the studied cases we have shown how the deterministic and the stochastic approaches can be used together for better understanding the system dynamics: the first gives information on the average time behaviour of the protocells population whereas the second elucidates the role of randomness.
In the case of homeostatic autopoietic vesicles, stochastic simulations have elucidated the role of random fluctuations that mainly affect the size distribution of aggregates in homeostatic conditions by selecting vesicles over an aggregation number threshold, since these vesicles are the most robust with respect to randomness of the surfactant degradation events.On the other hand, in the second case studied a deterministic general formula, the osmotic synchronization law, has been introduced in order to predict the occurrence of a protocell stationary self-production regime.It should be clear that this formula is strictly valid only for the assumptions inherent in our modelling approach and as such, this formula will be applicable to real-world vesicle systems only when the assumptions of our model are satisfied.
The more complex case of ribocells was then described.The ribocell is a hypothetical model of a true MAC since in principle it could exhibit the traits of Darwinian evolution.Here the attention focused on the capability of self-reproduction, i.e. the possibility of synchronization between membrane self-production and the genome self-replication.Deterministic analysis leads to the definition of optimal conditions for observing the occurrence of a stationary self-reproduction regime.Stochastic simulations performed in optimal conditions showed that the expected synchronization between the genome and membrane dynamics is highly affected by random fluctuations.As a consequence of this, the expected synchronization occurs only by chance and the behaviour of a Ribocells population is very far from that predicted deterministically, so that, generation by generation, a plethora of protocells-each one differing in the genetic material inherited from the mother-can be generated, and the initial ribocell population eventually dies through ribozymes segregation.
Are these models realistic?That is, do they have experimental counterparts?The first example (homeostatic autopoietic vesicles) was indeed realized in the laboratory by using simple chemical reactions to produce and degrade the vesicle constituents, which were fatty acids [12].The second study is instead purely theoretical and awaits possible experimental verification.The third model, which was in turn based on a hypothesis [36], also misses an in vitro model for a direct investigation of the time behaviour emerging from our simulations.
From these considerations comes the question on how to progress in vitro, towards (bio)chemical MAC models of increasing complexity.Several groups are currently involved in these pioneering attempts, following modern synthetic biology approaches for designing and assembling synthetic compartments.Previous work shows that the technology is already quite mature [8], also for the biochemical synthesis of lipids [48].The recent introduction of the so-called "droplet transfer" method is quite promising in this respect [49] (Carrara et al. manuscript in preparation), as well as its combination with microfluidic technologies [50].This method allows a facile encapsulation of several macromolecules inside a single compartment-a phenomenon that is generally hindered by the adverse statistics of co-encapsulation [51,52].Anchoring enzymes to a scaffold polymer [53] might also help the co-encapsulation of all required molecules to reconstitute biochemical paths (possibly taking advantage of enzyme proximity).Clearly, when more complex networks are reconstituted inside autopoietic or non-autopoietic vesicles, more sophisticated, and possibly stochastic biochemical models are required [54].
In conclusion, we would like to emphasize the common aspects of analysis and modelling of these (and other) systems, namely the need for a systemic approach that integrates (and couples) the internal reactions, the membrane dynamics, and the environment.This is perhaps the most important scientific message that emerges from numerical simulations of these complex systems.Since numerical modelling is carried out using true physical constants for all elementary molecular steps, it follows that simulation outcomes could provide quantitative data to guide experimentalists in the design and the construction of protocells or artificial cells for nano-technological applications.Finally, as demonstrated in this review, the use of stochastic modelling in conjunction with deterministic approaches can uniquely reveal intriguing dynamics in micro-compartmentalized complex multi-molecular systems and greatly helps to evaluate and understand basic mechanisms at the roots of biological behaviour.

Figure 1 .
Figure 1.Relationship between artificial chemical systems aimed at mimicking cellular behavior.Autopoietic systems are all the systems that fulfill the definition given by Maturana and Varela[7] and they can belong to completely different fields, like natural science, sociology or economy.A subset of this collection is represented by autopoietic chemical systems.This subset encompasses reverse[9] and direct[10] micelles (not explicitly shown) along with autopoietic vesicles[11].Depending on the kinetic regime, autopoietic vesicles can self-reproduce, stay in a homeostatic regime or decay[12].MACs can be classified as a subset of autopoietic self-reproducing vesicles since they may also exhibit the capability to evolve.On the other hand, we may call "self-producing" vesicles those vesicles that are able to self-produce at least the membrane, owing to an internal reaction (see the gray set).Whereas autopoietic vesicles are necessarily self-producing, self-producing vesicles are not necessarily autopoietic (their growth being limited to the membrane constituents).According to the Ruiz-Mirazo definition[13] all these chemical systems can be generally indicated as protocells.The three systems discussed in this article are explicitly shown in the diagram (in violet).

Figure 3 :
Figure 3: Protocell stability as a function of the reduced surface .

Figure 4 :
Figure 4: Time evolution of a reacting vesicle monitored by the grow control coefficient:(a)  > 3/2 osmotically-stressed growth, i.e. the volume increases faster, then it will reach an elastic tension condition and, above the limit of elasticity of the membrane, this will lead the vesicle to osmotic burst ( < 1−); (b)  = 3/2 continuous spherical growth, i.e. a spherical vesicle will increase its size without any change of shape ( = 1); (c)  < 3/2 reproductive growth, i.e. the surface increases faster than the two previous cases, the growing vesicle will become deflated, assuming a non-spherical shape (ellipsoidal, elongated or, generally speaking, a prolate shape) and the energy of the membrane will be higher due to a bending tension.

Figure 5 .
Figure 5. Autopoietic vesicles: Kinetic mechanism and different regimes on the left, schematic drawing on the right.The surfactant (lipid) precursor P is taken up by the vesicle, and it is converted to surfactant S with a generation rate v G .S can also react with an oxidant Y to give the byproduct W (with a degradation rate v D ).

Figure 6 .
Figure 6.Stochastic simulations of autopoietic vesicles in homeostatic conditions: (a) evolution of the vesicle size distribution (to each size class belong ennamers with size 2 m−1 < i ≤ 2 m except for the first class m = 1 where only monomers are included); (b) time evolution of a single aggregate made by an initial number of monomers equal to 1000.Reproduced with kind permission of the authors[35] and under the conditions of IOP Publishing.

Figure 7 .
Figure 7. Self-producing enzymatic vesicle in a growth regime.Only the surfactant (lipid) production takes place within the vesicle, whereas the catalyst (the enzyme) is not produced.P: surfactant precursor; S: Forming membrane surfactant; E: Enzyme, and W: Waste molecule.

Figure 8 .
Figure 8. Self-producing enzymatic vesicle deterministic curves (black lines and data) and stochastic simulation results (gray lines and data) comparison: time evolution of the core volume (left plot); division time against generations (right plot).Horizontal dashed lines represent values calculated using Equations (15) and(16).

Figure 10 :
Figure 10: Deterministic time behavior of the ribocell in optimal external conditions:[NTP] Ex = [P] Ex = 0.01M.At the starting time the genome was composed by 100 dimers of of R c R L and R c R P and the radius was 100 nm and the core volume 4.2  10 6 nm 3 .On the upper plots the time courses of the volume and the surface of the mother protocell are shown, while on the lower plot, the total number of lipase and polymerase strands present in the aqueous core as dimers and free monomers are reported against time.

Figure 11 :
Figure 11: Stochastic behavior of a population of 320 ribocells: population composition against time (on the left), schematic drawing of different protocells as result of vesicle division and random RNA strands distribution.Reproduced from[33] with kind permission of the author.Different protocells can emerges due to ribozymes segregation among the daughter protocells at the moment of mother protocell division: vesicles without any ribozymes (empty vesicles) or containing one lone R P or many filaments of c R P and/or c R L (inert vesicles), vesicles encapsulating R L (self-producing), vesicles containing more than one molecule of R P or both R P and c R P filaments (self-replicating genome), aggregates containing one R P filament and R L / c R L strands (reduced ribocells).

Table 1 :
Kinetic parameters for the in silico ribocell model at room temperature.