## Abstract

This paper is a theoretical attempt to gain insight into the problem of how self-assembling vesicles (closed bilayer structures) could progressively turn into minimal self-producing and self-reproducing cells, i.e. into interesting candidates for (proto)biological systems. With this aim, we make use of a recently developed object-oriented platform to carry out stochastic simulations of chemical reaction networks that take place in dynamic cellular compartments. We apply this new tool to study the behaviour of different minimal cell models, making realistic assumptions about the physico-chemical processes and conditions involved (e.g. thermodynamic equilibrium/non-equilibrium, variable volume-to-surface relationship, osmotic pressure, solute diffusion across the membrane due to concentration gradients, buffering effect). The new programming platform has been designed to analyse not only how a single protometabolic cell could maintain itself, grow or divide, but also how a collection of these cells could ‘evolve’ as a result of their mutual interactions in a common environment.

## 1. Introduction

To date, most theoretical approaches in the field of origins of life or artificial life have focused on the dynamics of ‘self-replicating’ entities, looking for the molecular roots of Darwinian evolution, regardless of the high complexity—and thus low prebiotic plausibility—of the real molecules, which would be involved in that kind of process. In contrast, the problem of how rather simpler molecular components and chemical processes could put together a protometabolic cellular system has received more occasional and marginal attention.

However, as the title of this special issue suggests, a new and ambitious research program on the artificial implementation of minimal—yet biologically relevant—cellular systems is beginning to take shape (e.g. Rasmussen *et al*. 2004; Szathmáry 2005). In this paper, we would like to contribute to that line of research by presenting a novel computational approach to model cell dynamics, which tries to bridge a gap between *in silico* and *in vitro* experiments, convinced that the final answer to the problem will come from the correct integration of both theoretical and experimental endeavours. Therefore, although we acknowledge the importance of ‘top-down’ approaches to the minimal cell project (e.g. Castellanos *et al*. 2004; Gabaldón *et al*. 2007), here we will address the problem only from a ‘bottom-up’ perspective, i.e. searching for the most simple components that could set up a chemical system with a potential to become a (proto)biological cell.

In order to do so, it is not enough to analyse the self-assembly processes that lead to vesicles or other types of supramolecular equilibrium structures, which form close, semi-permeable compartments (e.g. Mayer *et al*. 1997). It is very important to investigate the compartments whose ‘building blocks’ also take active part in chemical reaction networks or transport processes. The main reason behind this is that ‘the chemical logic of a minimum protocell’, as Morowitz *et al*. (1988) rightly claimed 18 years ago, must be the logic of a thermodynamically open system, which exchanges matter/energy with the environment so as to maintain a far-from-equilibrium dynamic organization.

Although one can find a good number of theoretical models that include the idea of cellular compartments (particularly relevant examples, from the field of origins of life, are Dyson 1982, 1999), they normally treat these as global constraints, which simply define the size of the system or the ‘limit number’ of components in it. A much more intricate problem is to investigate how the compartment itself may change its properties due to the processes taking place in its internal core—or in the surrounding environment—and how this can affect, in turn, the metabolic network. Unfortunately, theoretical or simulation models that deal with those aspects of protocellular systems are not so widespread in the literature.

Nevertheless, some relevant steps forward in that direction were already taken, for instance, with the work carried out on minimal ‘autopoietic’ systems, both at a theoretical level (Varela *et al*. 1974; McMullin & Varela 1997) and in connection with real experiments (Mavelli & Luisi 1996; Mavelli 2003). Computer simulations of Ganti's ‘chemoton’ model (Ganti 1975, 1987, 2004) can also be considered as contributions to a similar goal (Csendes 1984; Fernando & Di Paolo 2004; Munteanu & Solé 2006). In addition, the approach developed by Lancet and colleagues (Segré & Lancet 2000; Segré *et al*. 2001), although it does not specifically tackle the problem of integrating a cellular membrane and a protometabolic network, provides an interesting theoretical framework where the evolutionary potential of amphiphilic assemblies can be studied (the so-called ‘lipid world’ scenario). Finally, spatially explicit (discrete lattice) models have also been recently presented with the aim to analyse the emergence and self-maintenance of reproducing cellular structures (Ono & Ikegami 1999; Madina *et al*. 2003), or the protometabolic reaction–diffusion mechanisms that may lead to the cell division (Macía & Solé in press; see also Macía & Solé 2007).

In any case, none of the previous models have taken into account in a minimally realistic way the structural and dynamic properties of a cellular compartment, consisting of a bilayered lipid membrane and a water core, and focusing on the active role it plays in the constitution and maintenance of a protometabolic network. A membrane, especially if it is biologically inspired, should not be conceived as an abstract boundary surface that merely defines the volume of the cell, preventing the diffusion of its components. Even in the most simple protometabolism, the cellular membrane must be a functionally very active part of the system, since it channels the flow of matter–energy necessary for its continuous self-construction (Ruiz-Mirazo & Moreno 2004). Thus, the main aim of this work is to explore the dynamic properties of realistic membranes, showing in particular how these can get coupled with chemical reaction networks, and elucidating the critical role they must have played in the organization of prebiological systems.

In order to do this, we adopt a stochastic dynamic approach, making use of a computational platform designed to perform Monte Carlo simulations of complex chemically reacting systems. The role of stochasticity and fluctuations in the dynamics of biological cells is drawing the interest of an increasing number of researchers (McAdams & Arkin 1997; Arkin *et al*. 1998; Barkai & Leibler 1999; Rao *et al*. 2002), so that other similar simulation platforms are starting to be developed (Lu *et al*. 2004; Webb & White 2005). Within this general framework, we present here a new program based on the well-accepted Gillespie direct method (Gillispie 1976, 1977) that makes possible the stochastic simulation of minimal cell models where chemical reactions are coupled with self-assembly and diffusion processes. Moreover, since this program can also simulate the division or collapse of a cell, it allows studying of both the formation and the evolution of whole-cell populations.

The theoretical background of the stochastic approach and the main features of our software platform are discussed in §2, the basic assumptions to study (proto)cell dynamics are outlined in §3, and in §4 (figure 1) the different reaction schemes investigated are described. We first consider the stochastic dynamics of close bilayer membranes when only reversible self-assembly processes take place, in order to determine the general conditions for the stability of empty protocells (see scheme 0 in figure 1). Then, stable protocells are set in far-from-equilibrium conditions, introducing the irreversible production of fresh lipid molecules. Three different scenarios of increasing complexity will be analysed: external production of lipids (see scheme 1 in figure 1); synthesis of lipids taking place at the boundary, within the actual cell membrane (see scheme 2 in figure 1); and finally an internal protometabolic reaction cycle that produces new lipid molecules and waste compounds (see scheme 3 in figure 1).

The outcome of our simulations (explained in detail in §5) represents particular examples of cellular systems with a boundary that can not only grow/shrink and divide/burst but also change its composition and thus its structural/functional properties; a boundary that is closer to a real lipid bilayer, consisting of amphiphilic molecules with a specific head surface and volume, where transport/reaction processes are allowed to occur, and which is subject to physico-chemical constraints, such as osmotic pressure (caused by concentration unbalance and free water flow). Yet, the model described in §2 will not try to mimic nature in every single molecular detail, providing an accurate biophysical description of membrane properties and dynamics. It is just meant to inform us, at an abstract—though experimentally testable—theoretical level, about the way in which a cellular compartment could help to build up, from rather simple prebiotic conditions (e.g. a ‘monomer world’ scenario; Shapiro 2002), increasing levels of molecular and organizational complexity.

## 2. Material and methods

The basic programming platform used for our simulations is an object-oriented (C++) environment that has been developed from a previous code recently presented by Mavelli (2003) and Mavelli & Piotto (2006). The earlier program, called Reactor, was devoted to the exact simulation of the master equation associated with a homogeneous chemically reacting system according to the Monte Carlo procedure presented by Gillespie in the 1970s (Gillespie 1976, 1977). For the theoretical background and a detailed description of this method, the reader is referred to Mavelli (2003) or to the original works of Gillespie. In this paper, only a brief sketch of the stochastic approach to chemical kinetics is provided.

Given a general reaction mechanism occurring in a homogeneous well-stirred space domain,where *r*_{j,ρ} and *p*_{j,ρ} are the stoichiometry coefficients and A_{j} (*j*=1, 2, …, *N*) is the molecular species. The stochastic kinetic theory states that the probability of a reaction *ρ* to occur in the next infinitesimal time-interval [*t*, *t*+d*t*] can be expressed as follows:where *c*_{ρ} is the probability coefficient and *n*_{j} is the number of A_{j} molecules present at time *t* in the reactive domain. Therefore, for each reaction, a probability density function *w*_{ρ}(**n**) is defined, which depends only on the state of the system , i.e. on the molecular numbers of species involved in the reactive event *ρ* according to its stoichiometry.

On the other hand, the deterministic approach defines a reaction rate in terms of molar concentrations [A_{j}],for each elementary step of the kinetic mechanism. The main difference between the two approaches lies in describing the system in terms of continuous variables (i.e. the macroscopic concentrations) instead of discrete variables (i.e. number of molecules). If the ergodic hypothesis holds and, in the thermodynamic limit (large molecular populations), the deterministic and the stochastic approaches converge on average, then the probability coefficients can be expressed in terms of the deterministic rate kinetic constants,where *N*_{A} is Avogadro's number, *V* is the domain volume, and *m*_{ρ} is the *ρ*-reaction molecularity. In order to obtain the time course of all the species in the reacting systems, a set of ordinary differential equations has to be solved in the deterministic approach. This set gives the rate change of each molar concentration in terms of the reaction rates. On the other hand, if the stochastic approach is adopted, then one has to deal with the so-called master equation, a finite difference partial differential equation that gives the rate change of the state Markov probability density function of the reacting system. In this case, the analytical solution can be obtained only in few cases and consists in average time behaviour along with the standard deviations for all species populations.

Since going into the details of the stochastic kinetic theory is out of the scope of this paper, we will remark only that, even if for high numbers of reacting molecules the stochastic treatment gives the same result on average as the deterministic one, for small populations, random fluctuations can drive the system towards an unexpected evolution.

As the master equation is very hard to solve, Gillespie introduced an iterative Monte Carlo procedure to exactly simulate the stochastic time evolution of a reacting system. Taking the summation of all *w*_{ρ}(**n**),the jump density probability can be obtained, i.e. the probability that the system goes away from state **n** in the infinitesimal time-interval [*t*, *t*+d*t*]. Gillespie showed that the function *W*(**n**) is related to the density probability that a certain waiting-time interval Δ*τ* occurs between two consecutive reactive events as follows:(2.1)

Therefore, the stochastic evolution of a reacting system can be seen as a sequence of waiting-time intervals followed by instantaneous reactive events that take place according to the probability(2.2)where is the jump vector related to reaction *ρ*. As a consequence of this, the iterative Monte Carlo procedure at each step *i* has to draw two pseudorandom numbers, 0≤*g*_{1}, *g*_{2}≤1, the first to calculate the waiting time randomly distributed according to equation (2.1),and the second to select the new incoming event according to the probability (equation (2.2)),

It is important to remark that, in this method, the evolution time of the process (calculated by the sum of all waiting time Δ*τ*_{i}) is a real physical time, with a measure unit that depends on the probability coefficients or the kinetic constants used.

The new platform, called Environment, has been designed to deal with reacting systems that are no longer completely homogeneous. Although a detailed description of this platform and its potential applications is currently being elaborated (Mavelli & Ruiz-Mirazo in preparation), we sketch here its main features, just to illustrate the potentiality of the approach. Essentially, the new program decomposes the global system into a collection of different reactive domains, each of which is assumed to be homogeneous. Concentration gradients can be established across the boundary of two neighbouring domains and molecules be exchanged between them by means of diffusion processes. The density probability that a molecule A_{i} can cross the boundary surface between two neighbouring phases is calculated as , where *k*_{i} is the diffusion coefficient, [*A*_{i,1}] is the number of molecules of species A_{i} in phase 1, and *S*_{1,2} is the area of the interface surface. Thus, along with molecular reactions, a new type of event, diffusion processes, can now occur in the system, after a certain waiting time.

In this way, we put forward a strategy to overcome some of the limitations of the standard ‘well-stirred flow reactor paradigm’, opening the way to simulate chemical dynamic systems with an intrinsically more complex organization rather than just a population of reacting molecules freely in solution. A similar framework that models a cell system as a collection of different homogeneous reacting domains has recently been presented by Morgan *et al*. (2004). Nevertheless, these authors describe the cell time behaviour deterministically, by solving a set of ordinary differential equations, instead of simulating the probabilistic master equation. Moreover, the cellular membrane surface has no independent dynamic behaviour, but directly correlated with the core volume, since the lipid self-assembly process is not explicitly taken into account (nor the concentration of molecules in the external environment).

## 3. Cell model

As mentioned previously, our approach to model cell dynamics is based on the distinction of different homogeneous reactive domains, or phases, in the global system. The ‘environment’ represents the common aqueous phase, where simple molecules and cellular compartments are contained. Each cell, in turn, is decomposed into two different reactive domains: a hydrophobic lipid phase, ‘the bilayered cell membrane’, and an aqueous internal pool, ‘the cell core’.

The number of cells—and thus of reactive domains—is not fixed, but may increase if one of them divides. In fact, the program can follow the evolution of both ‘parent’ and ‘daughter’ cells after the division. The number of cells might also decrease if one of them bursts due to an osmotic crisis. In this case, the membrane opens up and all the internal water pool content is released to the global environment. Therefore, the cell immediately transforms into an open bilayer that behaves as a single hydrophobic domain and continues exchanging molecules with the water solution.

In each reactive domain, different types of molecular species and reaction processes are defined, according to the specificities of the system under investigation (see §4). Initial molecular concentrations and rate constants (i.e. reaction probability coefficients) are set as parameters in each domain before starting the simulations. Similarly, the permeability rules and value of the diffusion constants are established from the beginning, so as to define how the different domains are going to interact via diffusion processes. Moreover, quite importantly, free flow of water is assumed among the global environment and the core of the cell(s), in order to ensure the isotonic condition,that is, the global concentration of substrates inside and outside the membrane of each cell is constantly balanced. Therefore, throughout the simulation, at the end of every iteration, the core volume *V*_{core} of each cell is then rescaled in the following way:simulating an instantaneous flux of water to balance the osmotic pressure.

Therefore, while the volume of the environment remains fixed to its initial value (typically, *V*_{env.}=5.23×10^{−16} dm^{3}, three orders of magnitude bigger than the core volume of a 50 nm radius spherical cell), the volume of the cell core can change during the simulation. The cell membrane will also have its own dynamic behaviour, since it continuously exchanges lipids with both the internal and the external aqueous phases. At any time, the membrane surface *S*_{μ} of a certain cell can be calculated by the summationwhere *α*_{i} is the hydrophilic head area of all the surface active molecules located on the membrane and the factor 0.5 takes into account the double molecular layer. However, the core volume and the membrane surface must satisfy some geometrical constrains, so that the cell is actually viable. In fact, the conditions for the division or burst of a cell will depend on the relationship between the membrane surface and the core volume, within the following limits.

The actual membrane surface must be bigger than the theoretical spherical surface that corresponds to the actual core volume, otherwise the cell bursts.

The actual surface of the cell must be smaller than the theoretical surface that corresponds to two equal spheres of half the actual core volume, otherwise the cell divides into two statistically equivalent daughter cells.

These two limits establish the conditions for stability or viability of a cell in our model, i.e. the range of possible states in which it will not break or divide. In mathematical terms, defining the parameter *Φ* (‘viability coefficient’) as the ratio between the actual cell surface and the surface of a sphere with the same volume , the conditions for cellular stability become . Moreover, in order to take into account the lipid membrane elasticity and flexibility, the program also includes a tolerance parameter (*ϵ*) that slightly enlarges the stability range as follows:

This tolerance parameter was set at 0.1 in all reported simulations, so that the final relationship becomes 0.9≤*Φ*≤1.386. When *Φ*>1.386, the cell divides into two twin spherical daughters. We are aware that this also entails a crude simplification, because all processes of division should not lead to equal daughter cells. But it seems to be the first and easiest way to avoid some further suppositions (e.g. constant spherical shape, division when the initial size is doubled). Besides, it is important to stress that, in these conditions, cell growth could be observed without ending up in a division process or, alternatively, there could be a cell division without growth (just by membrane deformation). In any case, if the overall concentration of the internal species increases much higher in relation to the growth of the membrane, then this is bound to cause the breaking of the cell due to water inflow (‘Donnan effect’).

All simulations were carried out with bilayered membranes composed of a generic lipid L, whose molecular volume and head surface area are equal to 1.0 nm^{3} and 0.5 nm^{2}, respectively (taken as typical values for an amphiphilic molecule). Moreover, for the sake of standardizing initial conditions, these membranes were assumed to be initially perfect spheres, even if this condition is rather close to the low critical threshold of the system (1−*ϵ*)=0.9.

The probability for the uptake of a molecule A from an aqueous domain (the environment or the core) to the lipid bilayer was calculated as the product of a kinetic constant *k*_{Aμ} times the aqueous molecule concentration [A]_{env./core} multiplied by the membrane surface area *S*_{μ}, whereas the backward process probability (i.e. the release of a membrane component to an aqueous domain) was assumed proportional to the corresponding kinetic constant *k*_{A} times the number of molecules of that component in a bilayer. In the case of the exchange of lipids from and to the membrane, upon equating the forward and backward process probabilities *k*_{Lμ}[L]_{aq.}*S*_{μ}=*k*_{L}*n*_{L,μ}, the lipid aqueous equilibrium concentration can be obtained as follows:by approximating *S*_{μ}≈*α*_{L}*n*_{L,μ}/2 and remembering that the membrane is a bilayer. Since we typically set *k*_{Lμ}=1.0 M^{−1} *t*^{−1} and *k*_{L}=0.001 *t*^{−1}, the lipid equilibrium concentration is [L]_{eq.}=0.004 M, both in the cell core and in the external environment. Thus, all the concentrations have to be considered as molar concentrations. Furthermore, it should be mentioned that, since we have not defined the measurement for the unit of time, in principle, all simulations should be reported against time 10^{3}/*k*_{L}, although in the plots we always used the arbitrary unit notation (arb. unit).

## 4. Scenarios

In our bottom-up approach, different scenarios of increasing complexity are considered, in order to explore the dynamics of protocellular systems under the general constraints and assumptions previously discussed.

### (a) Scheme 0: equilibrium empty cell dynamics

Initially, we study the stochastic dynamics of empty cells (or bare closed bilayers), consisting of a single type of lipid, when only reversible self-assembling processes take place (see scheme 0 in figure 1). The aim is to determine the general conditions for the stability of closed membranes around equilibrium, so as to have a solid starting point for the next (far-from-equilibrium) simulations.

### (b) Scheme 1: irreversible lipid synthesis in the external environment

It is clear that, in order to obtain a more interesting (and biologically relevant) range of dynamic behaviours, the system must be progressively taken away from equilibrium conditions, introducing some process of production of lipids. Therefore, here we analyse the most elementary case, in which the membrane-reversible self-assembly is coupled with an irreversible reaction to produce lipids outside the cell (see scheme 1 in figure 1).

### (c) Scheme 2: minimal self-producing cells (the ‘Luisi’ scenario)

Now a separate pure organic phase is assumed to exchange a hydrophobic compound Z reversibly with the water environment, which is rapidly absorbed by the bilayer. Once in the membrane it spontaneously converts into L, promoting the autocatalytic growth and subsequent reproduction of the vesicle. In order to counterbalance this process, we also consider the possibility that the lipid is degraded into two simpler compounds Q, eventually released to the environment, which behaves as a sink for Q. These nonlinear autocatalytic conditions (see scheme 2 in figure 1) roughly correspond to those in which Luisi and colleagues performed their experiments with ‘self-(re-)producing’ vesicles (Walde *et al*. 1994). In fact, Luisi (1993) himself elaborated a simple model to account for a situation that also included the decay of the surfactant (as we do here), and later checked it experimentally (Zepik *et al*. 2001). This is why we refer to this general scheme as the Luisi scenario.

### (d) Scheme 3: ‘proto-chemoton’ cells

To advance our simulations and get closer to a protometabolic cell (although this would deviate us away from real experimental conditions), we considered that the chemoton model (Ganti 1975, 1987, 2004) could give a good insight, provided that we managed to adapt it, somehow, to the background conception from which this paper is written. In fact, although we regard the chemoton model as a very interesting scheme to study an intermediate stage in the process of the origin of life (e.g. a kind of ‘polynucleotide world’ embedded in a cellular protometabolism), the purpose of this work was to analyse previous transitions, which involved a lower level of molecular complexity (so modular ‘templates’, for instance, could not yet be there) and a higher level of noise and fluctuations (originally not present in the neat chemical design of Ganti).

Therefore, the final decision was to go for a less complicated reaction scheme, also inspired by Ganti's work (Ganti 1987, 2002): a sort of simplified version of the chemoton, which does not include the ‘template subsystem’. Therefore, this proto-chemoton consists of only two coupled autocatalytic subsystems: the membrane and the protometabolic network. The general idea is to carry out a computational analysis of that basic kinetic mechanism (scheme 3 in figure 1), putting it under realistic conditions (i.e. introducing membrane-reversible self-assembly, diffusion processes, osmotic pressure, fluctuations).

## 5. Results and discussion

### (a) R-scheme 0: equilibrium empty cell dynamics

The first relevant results of our simulations concern the role of random fluctuations as a possible factor of instability for a cellular system. We mimic the case of a vesicle suspension in pure water solution or in the presence of an osmotic buffer, which is a compound that cannot cross the membrane but has the same concentration inside and outside the cell, for instance an inorganic salt. As shown in figure 2, spherical vesicles in equilibrium with an aqueous concentration of lipid can fluctuate around the initial value of the volume, but below a critical value of 30 nm, vesicles underwent an osmotic burst. Figure 3*a* shows the lipid concentration in the internal aqueous solution, which fluctuates around the expected equilibrium value of 0.004 M. In figure 3*b*, the trend of *Φ* for an unstable cell is reported. The influence that the presence of an osmotic buffer has on vesicle stability was then investigated. Moreover, the results portrayed in figure 4 clearly show that as the buffer concentration increases, the oscillations in the volume remarkably decrease, since the osmotic pressure across the membrane is less unbalanced by stochastic fluctuations of the lipid concentration in the core.

We also checked the response of a stable vesicle to an external perturbation: the addition of some osmotic buffer molecules or fresh lipids from outside. When the buffer is externally added to pre-existing vesicles in a pure water solution, a kind of ‘osmotic shock’ occurs (data not shown). Indeed, if the final concentration of the added buffer [B]_{env.}≫[L]_{eq.} and [B]_{core}=0, then this produces an immediate shrinkage of the cell (division followed by osmotic crisis). But if the addition is done to an already buffered solution [B]_{env.}>[B]_{core}≫[L]_{eq.}, then a gentle deflating of the core volume can be induced and the vesicles remain stable but no more spherical: .

The external addition of fresh lipids is a more interesting case, since they can be integrated in the membrane, provoking growth and division of the initial cell. For these simulation runs, a relatively high concentration of buffer was set ([B]_{env.}=[B]_{core}=0.02 M) in order to prevent the previously described osmotic shock.

In figure 5*a*,*b*, we show a typical result obtained by adding to a 50 nm radius cell (fluctuating around equilibrium) an amount of lipid 10 times higher than the amount present in the initial cell, which corresponds to a concentration increment Δ[L]=0.004 M. After the addition (indicated by the vertical arrow in the graph), an instantaneous decrease in the core volume due to the osmotic shock (figure 5*a*, black curve) is observed, followed by an increase in the surface (figure 5*a*, grey curve) at a constant volume value. As soon as *Φ* reaches the upper limit for stability, the original cell splits (figure 5*b*, grey curve), and this happens twice again during the simulated time-interval. Looking at the time evolution of the total number of cells (figure 5*b*, black curve), one can note that the first division step is not followed by an increase in the cell number, and infer that the daughter cell must have ‘died’ very soon due to an osmotic crisis. Instead, the second and third divisions are normal which take place with an increase in the population (from one to two and two to four), but shrinking to smaller cells that are less stable and quickly disappear. In contrast, when the amount of added lipids per cell is smaller and the buffer concentration is higher, the growth is slower but the percentage of successful divisions clearly increases. As an example, figure 5*c* reports the simulation carried out starting with a population of 10 identical cells in high buffer concentration ([B]_{env.}=[B]_{core}=0.2 M) and adding the same amount of lipid as before. This leads exactly to double the cell population. In conclusion, the general tendency when fresh lipids are added from the outside is an increase in the number of cells but with a mean size reduction (progressive shrinking).

### (b) R-scheme 1: irreversible lipid synthesis in the external environment

In order to overcome the previous setting and investigate non-equilibrium cell dynamics, the reversible self-assembly processes have to be coupled with irreversible chemical reactions. The simplest case is to add to the external environment a lipid precursor X that can spontaneously turn into lipid. In this way, the release of extra L molecules to the environment can be tuned by the value of the kinetic constant associated with the irreversible process *k*_{XL} [*t*^{−1}] (figure 1). When *k*_{XL} is very high (compared with the self-assembly kinetic constants: *k*_{XL}≫*k*_{Lμ}[L]_{eq}≫*k*_{L}), the same time behaviour as in the case of an instantaneous addition of lipids is expected. Therefore, simulations were carried out in the same conditions of the lipid addition simulations (previous case), starting with a 50 nm radius cell and setting [X]_{0}=0.004 M, equal to the concentration increment Δ[L]. As shown in figure 6, when the kinetic constant is high (*k*_{XL}=10^{6}×*k*_{L}), the behaviour converges to the ‘limit case’ and this is a good check for the coherence of our program. But if *k*_{XL} is sufficiently low (*k*_{XL}=10*k*_{L} or below), the growth is also slower, as expected. In any case, if the simulation is left to run long enough, subsequent divisions with progressive shrinking of the population mean size could not be avoided. This is due to the fact that there are no reasons for an increase in the core volume, so the surface growth leads to cell splitting with size reduction.

### (c) R-scheme 2: minimal ‘self-producing’ cells

In the case of scheme 2 in figure 1, we first considered the situation without the degradation or decay of the lipid. The typical conditions under which simulations were carried out are those of an initially spherical vesicle of 50 nm in equilibrium with lipid molecules ([L]_{env.}=[L]_{core}=0.004 M and *k*_{L}=0.001 *t*^{−1}, *k*_{Lμ}=1.0 M^{−1} *t*^{−1}) at a relatively high buffer concentration to ensure stability ([B]_{env.}=[B]_{core}=0.2 M). The density probability of the Z release from the γ pure organic phase into the water environment is defined as *k*_{Z}*S*_{γ}[Z_{γ}], where [Z_{γ}] is the phase molar density and *S*_{γ} is the interface area between the organic and the water environment. The γ-phase is considered as a macroscopic reservoir of the hydrophobic compound Z that floats on top of the water phase, as described by Luisi and co-workers (Walde *et al*. 1994). Thus, the release density probability is assumed to be constant and set equal to 10^{−4} *t*^{−1}, whereas the density probability of the Z phase separation from the water environment is calculated by *k*_{Zγ}*S*_{γ}[Z]_{out}. We establish *k*_{Zγ}*S*_{γ}=1.0 M^{−1} *t*^{−1}, so that the equilibrium solubilization of Z in the water environment is [Z]_{eq.}=10^{−4} M, in agreement with the assumption that Z is water insoluble. In addition, since the density probability of the Z uptake from the aqueous environment to the cell membrane is *k*_{Zμ}*S*_{μ}[Z]_{env.}, we assume *k*_{Zμ}=10^{4} to mimic a highly spontaneous process. Fixing these parameter values as mentioned, a set of different simulations was run changing *k*_{ZL}, the kinetic constant of the irreversible transformation Z_{μ}→L_{μ} (which takes place within the cell membrane). Figure 7 shows a typical time behaviour of the system in one of these runs, obtained for *k*_{ZL}=0.001 *t*^{−1}. What is usually observed is a decrease in the average volume and membrane surface of the cells, as a consequence of the successive divisions, while their population grows. Increasing *k*_{ZL} from 0.001 to 0.1 does not alter the overall results, but increases only the division frequency. It should also be noted that the volume fluctuations are practically suppressed by the high buffer concentration (at the scale used in figure 7*a*).

Later, scheme 2 including the irreversible lipid breakdown into two equal molecular fragments Q (then released to the water environment) was analysed. Setting *k*_{Q}=*k*_{L}=1.0, three altogether different outcomes were observed as a function of the *k*_{LQ} associated with the process L_{μ}→2Q_{μ}. In figure 8*a*, the time courses of *Φ* are reported in these three cases: (i) rapid osmotic crisis (*k*_{LQ}=3.0×10^{−6}), (ii) homeostatic regime (*k*_{LQ}=2.0×10^{−6}) and (iii) continuous cell division with an increase in the population (*k*_{LQ}=1.0×10^{−6}). Figure 8*b* shows the particular case of the homeostatic regime, where the surface and the volume of the cell reach a stationary state and fluctuate around it. It should be noted that the cell membrane is not a perfect sphere and it happens to be in tension since 0.9<*Φ*<1.0. In figure 8*c*, the case with subsequent division steps is portrayed, showing that here also the increase in population corresponds to a shrinkage of the average cell size.

### (d) R-scheme 3: proto-chemoton cells

The final scenario explored with our simulation program is the case of a cell that is able to produce lipids, thanks to the internal protometabolic cycle (see scheme 3 in figure 1). Introducing Ganti's (2002) simplified scheme in our way to model cell dynamics involves the following boundary conditions.

Regarding permeability rules, the precursor molecule X and waste product W are allowed to cross the membrane (depending on their particular concentration gradient at each time-step). Lipid molecules could also cross from the core to the environment (or vice versa) but in two steps: joining the membrane first and then leaving it to go to the aqueous solution on the other side. Lipids L—together with some buffer compound B—will be present both inside and outside the cell. Instead, the components of the metabolic cycle A

_{i}are assumed to be completely impermeable substances (only present in the internal core).Regarding concentrations, the environment is supposed to be big enough, so that it works as an infinite source of precursors ([X]

_{env.}=const.) and a sink for waste products ([W]_{env.}=const.). As for the initial values of a typical run, [L]_{env.}=0.004 M=[L]_{core}(equilibrium concentration), [B]_{env.}=[B]_{core}=0.2 M, [X]_{env.}=0.001 M, [X]_{core}=0 M, [W]_{env.}=0 M, [W]_{core}=0 M, [A_{1}]_{core}=0.001 M and [A_{i}]_{core}=0 M, where*i*≠1.Regarding the kinetic constants, we consider

*k*_{i}=1.0 (for all the forward reactions in the cycle) and (for the backward ones). This asymmetry amounts to say (Csendes 1984) that the X component (i.e. the ‘food’ or ‘energy input’ into the cycle) is a high free energy compound, compared with all the other compounds of the cycle (including the W and L ‘by-products’). In other words, the system is somehow forced to run in non-equilibrium conditions, so the metabolic cycle, in practice, is not completely reversible.

Under these conditions, we studied the time behaviour of the cellular system as a function of the initial size and the diffusion constants *k*_{W} and *k*_{X}. We found out (as portrayed in figure 9*a*) that, in this scenario, there is also a critical size below which the cell undergoes an osmotic crisis and bursts. For a fixed value of the kinetic constants (*k*_{W}=0.1 and *k*_{X}=0.01), vesicles with *R*≤40 nm were found unstable. Moreover, above that minimal size threshold, the bigger the cell, the faster it grows and divides. In figure 10, a subsequent division regime is reported for an initial 50 nm radius cell (with the same *k*_{W}=0.1 and *k*_{X}=0.01). The number of cells formed at the end of this simulation was 16 and no osmotic crisis occurred. In figure 10*b*, the core concentrations of the parent cell are shown. But it is perhaps more interesting to note, from figure 10*a*, that this is the first case in our *in silico* experiments in which the increase in cell population corresponds to an increase in the average cell size. This is due to the metabolic cycle that increases the overall core concentration, forcing water into the cell in order to fulfil the isotonic condition.

Finally, figure 9*b* shows another important finding in this context: given an initial stable cell (e.g. a 50 nm radius spherical cell), the critical parameter that determines whether it is going to divide or burst by osmotic crisis is *k*_{W} (i.e. the diffusion constant of the waste product). Irrespective of the value of *k*_{X} (that seems only to put forward/delay slightly the entry of X to the compartment, and therefore the time for the first division), for *k*_{W}=0.001 or below the system underwent a crisis, whereas for higher values it grew and divided. Hence, we can conclude that the viability of the cell strongly depends on its capacity to get rid of its waste material.

## 6. Final remarks and outlook

Under the general assumptions of our cell modelling approach and the specific reaction networks explored, we found particular parameter values that led to the emergence of biologically interesting time behaviours, including subsequent reproductive cycles of the initial cells. The different simulation systems and settings analysed here were meant to be close to real experimental conditions and, indeed, produced congruent results. We also explored a more complex and abstract situation (the case of proto-chemoton cells) that points in the direction where, from our perspective, future research should be conducted. It was quite significant that only in this final case, when the lipid production takes place within the cell core, we observed reproduction with actual growth in the size of the cellular system. Nevertheless, there were other remarkable results, like the homeostatic case obtained in the Luisi scenario, somehow analogous to the results obtained in real *in vitro* experiments (Zepik *et al*. 2001).

However, we must point out that the present paper constitutes just a first step in a hopefully long research avenue opened with this new approach to the simulation of protocell dynamics. On the one hand, although a relatively wide variety of behaviours were already observed, we still expect a richer dynamic landscape for the same—or roughly similar—minimal cell scenarios. In fact, a more thorough and statistically well-grounded analysis of the parameter space for the different schemes introduced is still required, specially with the aim to get closer to critical conditions (e.g. very low populations of metabolites), where the role of fluctuations and stochasticity will become more relevant.

On the other hand, this computational platform has many potential applications to study other minimal, infrabiological cellular systems, as well as some aspects of fully-fledged living cells. At present, we consider that the most natural continuation of the work carried out so far would be an elaboration of the program, so that it can include simple oligomers, like peptide chains, in the scheme. This would allow us to explore the role that such new components (with their own molecular characteristics, autocatalytic dynamics, etc.) could play in stabilizing the protocellular compartments, forming rudimentary transmembrane channels that improve the control of osmotic unbalance or the capacity to throw away waste products. In this sense, we search for minimal lipid–peptide cell models that could help us understand the origins of more autonomous molecular agents (Kauffman 2003; Ruiz-Mirazo & Moreno 2004).

In a similar way, the introduction of new programming objects in the platform to account for reaction schemes with more complex types of molecule (like oligonucleotides or even RNA fragments) would lead to further simulation work that could be helpful to tackle the problem of integrating protometabolic cells with molecular replication mechanisms. In this new scenario, a very important aspect to examine would be how these replication mechanisms can contribute to a more reliable reproduction of the whole cellular system. Analogously to what we did here with the simplified proto-chemoton scheme, the complete chemoton model of Ganti (1975, 2004; i.e. including the template autocatalytic subsystem) could then be tried under realistic cell conditions.

Thus, there is ground for a number of very promising developments of the present platform. We just hope that this and our possible future work together will illuminate some of the critical issues researchers will have to face on their way to achieve the implementation of artificial autonomous cells.

## Acknowledgments

K.R.-M. carried out this research, thanks to a fellowship from the Basque Government and the research grant 9/UPV 00003.230-15840/2004. The authors would also like to acknowledge the financial help from COST (Action D27) and the HPC-Europa program, which made possible their collaboration.

## Footnotes

One contribution of 13 to a Theme Issue ‘Towards the artificial cell’.

- © 2007 The Royal Society