Royal Society Publishing

Stochastic simulations of minimal self-reproducing cellular systems

Fabio Mavelli , Kepa Ruiz-Mirazo


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).

Figure 1

Scheme 0: equilibrium empty cell dynamics. Scheme 1: irreversible lipid synthesis in the external environment. Scheme 2: minimal ‘self-producing’ cells. Scheme 3: ‘proto-chemoton’ cells.

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,Embedded Imagewhere rj,ρ and pj,ρ are the stoichiometry coefficients and Aj (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+dt] can be expressed as follows:Embedded Imagewhere cρ is the probability coefficient and nj is the number of Aj 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 Embedded Image, 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 [Aj],Embedded Imagefor 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,Embedded Imagewhere NA 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),Embedded Imagethe 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+dt]. 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:Embedded Image(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 probabilityEmbedded Image(2.2)where Embedded Image 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≤g1, g2≤1, the first to calculate the waiting time randomly distributed according to equation (2.1),Embedded Imageand the second to select the new incoming event according to the probability (equation (2.2)),Embedded Image

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 Ai can cross the boundary surface between two neighbouring phases Embedded Image is calculated as Embedded Image, where ki is the diffusion coefficient, [Ai,1] is the number of molecules of species Ai in phase 1, and S1,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,Embedded Imagethat 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 Vcore of each cell is then rescaled in the following way:Embedded Imagesimulating an instantaneous flux of water to balance the osmotic pressure.

Therefore, while the volume of the environment remains fixed to its initial value (typically, Venv.=5.23×10−16 dm3, 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 summationEmbedded Imagewhere α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.

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

  2. 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 Embedded Image, the conditions for cellular stability become Embedded Image. 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:Embedded Image

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 nm3 and 0.5 nm2, 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 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 kA 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]aq.Sμ=kLnL,μ, the lipid aqueous equilibrium concentration can be obtained as follows:Embedded Imageby approximating SμαLnL,μ/2 and remembering that the membrane is a bilayer. Since we typically set k=1.0 M−1t−1 and kL=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 103/kL, 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 3a shows the lipid concentration in the internal aqueous solution, which fluctuates around the expected equilibrium value of 0.004 M. In figure 3b, 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.

Figure 2

Stability of spherical membranes with different radius (R) in a pure water solution (cell volume versus time): the internal volume fluctuates continuously around the geometrical value 4/3πR3. Fluctuations can bring smaller structures (R≤30) to collapse due to an osmotic crisis.

Figure 3

(a) Lipid concentration in the intracellular water pool reported against time in the case of an initial spherical membrane of 50 nm radius: [L]core fluctuates around the equilibrium value 0.004 M. (b) Φ factor (Φ=S/Ssphere) reported against time for an initial spherical membrane of 30 nm radius. Φ fluctuates around 1.0 (i.e. the value of a perfect sphere). An osmotic crisis occurs when the dashed line is crossed; that is, when the deviation from the threshold value goes beyond 10% (tolerance: ϵ=0.1).

Figure 4

Effect of the osmotic buffer on the fluctuations of core volume. As buffer concentration increases, fluctuations decrease. Osmotic buffer effect on 25 nm size membrane: (a) [B]out=[B]core=0.05 M and (b) [B]out=[B]core=0.5 M. (c) Osmotic buffer effect on membrane stability: the average volume fluctuations are reported against the buffer concentration.

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: Embedded Image.

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 5a,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 5a, black curve) is observed, followed by an increase in the surface (figure 5a, grey curve) at a constant volume value. As soon as Φ reaches the upper limit for stability, the original cell splits (figure 5b, grey curve), and this happens twice again during the simulated time-interval. Looking at the time evolution of the total number of cells (figure 5b, 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 5c 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).

Figure 5

Addition of fresh lipids to equilibrium cells. (a,b) A single 50 nm radius cell at medium buffer concentration: [B]env.=[B]core=0.02 M. The amount of added lipids was 10 times of that required for the initial cell: [L]Tot=4.4×10−3 M before the addition and [L]Tot=8.4×10−3 M after the addition. (c) The addition of the same amount of lipids as before but to a suspension of 10 identical 50 nm cells at higher buffer concentration [B]env.=[B]core=0.2 M.

(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 kXL [t−1] (figure 1). When kXL is very high (compared with the self-assembly kinetic constants: kXLk[L]eqkL), 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 (kXL=106×kL), the behaviour converges to the ‘limit case’ and this is a good check for the coherence of our program. But if kXL is sufficiently low (kXL=10kL 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.

Figure 6

Irreversible lipid synthesis in the water environment: different time courses of the surface membrane for different values of the kinetic constant kXL (see scheme 1 in figure 1). For low values of kXL, the growth of the system is slow (as expected), while if kXL increases the time behaviour tends to the limit case of the instantaneous lipid addition. Simulations were carried out starting from a single 50 nm radius cell in equilibrium with lipids at relatively high buffer concentration [B]env.=[B]core=0.2 M.

(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 kL=0.001 t−1, k=1.0 M−1t−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 kZSγ[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−4t−1, whereas the density probability of the Z phase separation from the water environment is calculated by kSγ[Z]out. We establish kSγ=1.0 M−1t−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 kSμ[Z]env., we assume k=104 to mimic a highly spontaneous process. Fixing these parameter values as mentioned, a set of different simulations was run changing kZL, 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 kZL=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 kZL 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 7a).

Figure 7

Minimal ‘self-producing’ cell without lipid decay (see scheme 2 in figure 1): (a) time course of the average core volume and average membrane surface of the cells and (b) time course of the average Φ value (Φ=S/Ssphere) and the cell population.

Later, scheme 2 including the irreversible lipid breakdown into two equal molecular fragments Q (then released to the water environment) was analysed. Setting kQ=kL=1.0, three altogether different outcomes were observed as a function of the kLQ associated with the process Lμ→2Qμ. In figure 8a, the time courses of Φ are reported in these three cases: (i) rapid osmotic crisis (kLQ=3.0×10−6), (ii) homeostatic regime (kLQ=2.0×10−6) and (iii) continuous cell division with an increase in the population (kLQ=1.0×10−6). Figure 8b 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 8c, 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.

Figure 8

Minimal ‘self-producing’ cell including the lipid decay (see scheme 2 in figure 1). (a) On the right axis, the average Φ versus time is reported for different values of the kinetic constant kLQ: 3.0×10−3kL (osmotic crisis), 2.0×10−3kL (homeostatic regime), and 1.0×10−3kL (division with shrinking); on the left axis, the increasing number of cells in the population versus time is reported (in this last case with division). (b) Surface and volume of the cell versus time in the case of the homeostatic regime. (c) Surface and volume of the cell versus time in the case of multiple division with shrinking.

(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.

  1. 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 Ai are assumed to be completely impermeable substances (only present in the internal core).

  2. 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, [A1]core=0.001 M and [Ai]core=0 M, where i≠1.

  3. Regarding the kinetic constants, we consider ki=1.0 (for all the forward reactions in the cycle) and Embedded Image (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 kW and kX. We found out (as portrayed in figure 9a) 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 (kW=0.1 and kX=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 kW=0.1 and kX=0.01). The number of cells formed at the end of this simulation was 16 and no osmotic crisis occurred. In figure 10b, the core concentrations of the parent cell are shown. But it is perhaps more interesting to note, from figure 10a, 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.

Figure 9

(a) ‘Proto-chemoton’ vesicles of different sizes (60–30 nm): time evolution results in terms of the Φ value (Φ=S/Ssphere). A critical value of R>40 nm was found to overcome an eventual osmotic crisis. The bigger the size, the higher will be the stability and growth rate (shorter time to the first division step). (b) Proto-chemoton vesicles of 50 nm size, for which the diffusion constant of the waste product (kW) was varied. It is evident from the graph that this is a fundamental parameter to guarantee the stability of the cell.

Figure 10

Typical results for 50 nm ‘proto-chemoton’ vesicles (kW=0.1 and kX=0.01). (a) The average cell volume and membrane surface are plotted against time. (b) The time evolution of the concentrations of all the internal species (of the original or ‘mother’ cell) is reported.

Finally, figure 9b 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 kW (i.e. the diffusion constant of the waste product). Irrespective of the value of kX (that seems only to put forward/delay slightly the entry of X to the compartment, and therefore the time for the first division), for kW=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.


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.


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


View Abstract