## Abstract

The reproduction of a living cell requires a repeatable set of chemical events to be properly coordinated. Such events define a replication cycle, coupling the growth and shape change of the cell membrane with internal metabolic reactions. Although the logic of such process is determined by potentially simple physico-chemical laws, modelling of a full, self-maintained cell cycle is not trivial. Here we present a novel approach to the problem that makes use of so-called symmetry breaking instabilities as the engine of cell growth and division. It is shown that the process occurs as a consequence of the breaking of spatial symmetry and provides a reliable mechanism of vesicle growth and reproduction. Our model opens the possibility of a synthetic protocell lacking information but displaying self-reproduction under a very simple set of chemical reactions.

## 1. Introduction

In 1952, Alan Turing published a very influential paper in *Philosophical Transactions of the Royal Society*, entitled ‘The Chemical Basis of Morphogenesis’. In that article, Turing proposed a possible solution to the problem of how developing systems can become heterogeneous, spatially organized entities starting from an initially homogeneous state (Turing 1952; Meinhardt 1982; Murray 1989; Lengyel & Epstein 1992). Turing showed that an appropriate compromise between local reactions and long-range communication through diffusion could generate macroscopic spatial structures. The interplay of both components was described in terms of a system of partial differential equations, so-called reaction-diffusion (RD) equations, namely a set(1.1)(1.2)Here, *C*_{1} and *C*_{2} indicate the concentrations of the two morphogens and their specific molecular interaction are described by the reaction terms *Φ*_{i}(*C*_{1}, *C*_{2}) with *i*=1, 2. These reactions could be, for example, activations, inhibitions or autocatalysis and degradation. The concentrations are spatial- and time-dependent functions, i.e. *C*_{1}=*C*_{1}(*r*, *t*) and *C*_{2}=*C*_{2}(*r*, *t*). Here *r* indicates the coordinates of a point *r*∈*Γ* where *Γ* is the spatial domain where reactions occur.

The last terms on the right-hand side of both equations stand for diffusion over space: here *D*_{1} and *D*_{2} are the corresponding diffusion rates, indicating how fast each molecule diffuses through space. If *Γ* is a two-dimensional area, we would have *r*=(*x*, *y*) and the diffusion term reads(1.3)which can be properly discretized using standard numerical methods. The key idea of Turing instabilities is that small initial fluctuations can be amplified through the reaction terms (typically nonlinear) and their effects propagate through space thanks to diffusion. These patterns are generated from an initially almost homogeneous distribution of morphogens. Specifically, we use the equilibrium concentrations as reference state obtained from the condition d*C*_{1}/d*t*=d*C*_{2}/d*t*=0. It is not difficult to see that if the system starts exactly from this homogeneous state, it will remain there forever. Now consider a very small perturbation of such homogeneous state, where the initial concentration is now with *ξ*(*r*, 0) a small noise term. Such initial fluctuations (inevitable due to the intrinsic noise) can be amplified by reactions and propagated by diffusion. Under some conditions (Turing 1952; Murray 1989), they can generate large-scale spatial structures with a characteristic scale. Such scale (wavelength) only depends on the intrinsic parameters involved in the RD terms. Such process of amplification of fluctuations can eventually shift the spatial distribution of morphogens from homogeneous to heterogeneous. Technically, this corresponds to a so-called symmetry breaking (SB) phenomenon: the initial symmetry defined by the nearly homogeneous spatial state is broken. The consequence of such SB is a heterogeneous pattern of morphogen concentrations (Nicolis & Prigogine 1977; Nicolis 1995).

Although the mechanisms underlying pattern formation in multicellular systems are typically richer than the previous RD scheme, they provide the appropriate framework to explain different situations. Some examples are pattern formation in fish (Kondo & Asai 1995), bacterial growth in two-dimensional cell cultures (Golding *et al*. 1998), sea shell patterns (Meinhardt 1998), the skin of vertebrates (Maini 2003; Suzuki *et al*. 2003), the self-organization of ant cemeteries (Theraulaz *et al*. 2002) and the spatial distribution of population densities in ecosystems (Solé & Bascompte 2006).

One particular scenario where living systems develop a spatial asymmetry is provided by single cells in morphogenesis. During early morphogenesis, cells often display a spatially asymmetric distribution of some molecules that appear to be preferentially located in different cell poles (Alberts *et al*. 2002). Such changes involve complex networks of molecular interactions and the reorganization of the cytoskeleton and are typically triggered by fertilization. At a simpler level, dynamical instabilities generating waves have been found in the cell cycle division of some bacteria. This seems to be the case in *Escherichia coli*, where a wave of protein concentrations, with rapid oscillations between the two membrane poles, seems to organize the division process (Raskin & de Boer 1999; Hale *et al*. 2001). Although the full mechanism is rather complex and involves polymerization processes beneath the cytoplasmic membrane, the mechanism driving the cycle is simple.

The problem of how supramolecular assemblies self-reproduce is at the heart of current efforts directed towards building synthetic protocells (Szostack *et al*. 2001; Rasmussen *et al*. 2005). Although many works have been devoted to studying the chemical coupling between vesicles, enzymes and information molecules, they have so far failed to produce reliable self-replicating protocells, although significant steps have been performed (Luisi *et al*. 1999; Oberholzer & Luisi 2002; Nomura *et al*. 2003; Takakura *et al*. 2003; Hanczyc & Szostak 2004; Noireaux & Libchaber 2004; Deamer 2005; Noireaux *et al*. 2005). One of the most promising approaches involves a top-down approximation using microscopic lipid vesicles incorporating pre-existing biological molecules (Noireaux & Libchaber 2004). The goal of this approach is finding a coupled set of reactions linking enzymes and/or information molecules with a container in such a way that growth and eventually reproduction can be achieved.

One problem with this top-down approach is that reliable processes leading to cell division need the appropriate coupling between the molecules involved (Solé *et al*. 2006). Such coupling must be able to increase cell size until some instability triggers vesicle splitting by generating a spatial breaking of symmetry through some active growth and deformation process. Properly designed, such a process can eventually end up in vesicle splitting.

Here, we propose a well-defined mesoscopic physico-chemical scenario of SB instabilities that is shown to generate the appropriate, self-maintained spatial heterogeneities.

## 2. Protocell model

The analysis of minimal cellular structures can contribute to a better understanding of possible prebiotic scenarios in which cellular life could have originated (Maynard Smith & Szathmáry 2001) as well as to the design and synthesis of artificial protocells (Rasmussen *et al*. 2004). Considering a minimal cell structure, basically two mechanisms can give origin to reproduction: spontaneous division and induced division. Spontaneous division takes place when the vesicle grows until splitting into two daughter cells becomes energetically favourable. Induced division is a more complex mechanism that allows for internally controlling the division process. These two scenarios are very different. In most models of protocell replication, reactions are well described but the container appears only implicitly defined (Kaneko & Yomo 2002; Gánti 2003; Kaneko 2005, 2006; Munteanu & Solé 2006; Sato & Kaneko 2006) and thus an important ingredient of protocell dynamics (the explicit presence of a changing container) is missing. So far, models dealing with some type of self-reproducing spatial structure have been limited to special type lattice systems (Ono & Ikegami 1999; Madina *et al*. 2003).

In an early work, the Russian biomathematician Nicolas Rashevsky (Rashevsky 1960; Solé *et al*. 2006) explored the conditions for instability-induced cell division. He concluded that during the process of membrane growth there is a critical radius beyond which spontaneous division is energetically favourable. Moreover, he suggested that time and space-variable osmotic pressures were one of the most suitable mechanisms inducing membrane division.

In a previous work, we have shown that osmotic-induced division is a feasible mechanism of vesicle self-replication (Macía & Solé submitted). In this framework, the non-uniform distribution of osmotic pressures along the membrane is related to the non-uniform, enzyme-driven metabolite distribution inside the vesicle, with metabolic reactions taking place in specific locations, where metabolic centres are located. These centres (using the term coined by Rashevsky) could be specifically designed, synthetic trans-membrane proteins. This method, however, can trigger just a single vesicle division cycle. After division, only one metabolic centre is present at each daughter cell and the division process cannot start again.

In our model (see below), replications take place indefinitely (provided that the appropriate precursors are available) and no enzymatic centres are required. A first approximation to such a minimal cell structure considers a continuous closed membrane involving some simple, internal metabolism. Here we propose a chemical mechanism coupled with vesicle growth, which generates the appropriate osmotic pressure distribution along the membrane (figure 1). These pressure changes are generated by the interaction of Turing-like instabilities with vesicle dynamics. They are able to induce the correct vesicle deformation and eventually cell division.

Although previous work on pattern formation has already considered changing spatial domains due to tissue or organism growth (Meinhardt 1982; Varea *et al*. 1997; Painter *et al*. 1999) and even membrane-bound Turing patterns (Levine & Rappel 2006), as far as we know, this is the first model where the boundaries are themselves a function of the RD dynamics, including membrane permeability and osmotic pressures altogether. As such, our model actually defines a totally new class of spatially extended dynamical system (figure 2).

## 3. Metabolism–membrane system

The problem of vesicle growth and division is not primarily thermodynamic, but kinetic (Morowitz 1992). If a reliable cycle of cell growth and splitting is to be sustained, we need: (i) precursors provided from the external environment and (ii) a restricted microenvironment where an appropriate set of reactions can drive the system out from equilibrium. Low-energy molecules and membrane precursors would be selectively transported across the membrane and high-energy compounds would be produced through energy-conversion processes.

Here, we explicitly define all the components of our protocell model. The main goal is to introduce a set of reactions to be represented by a set of *n* RD equations, namely(3.1)with *i*=1, …, *n* the index associated to the *i*th morphogen. However, we will extend the formalism by incorporating a changing boundary that now acts as a permeable membrane, also coupled to the reactions described by *Φ*_{i}. These reactions will define the protocell metabolism. Since osmotic pressures are associated with differences in molecular concentrations, active mechanisms generating spatial heterogeneity are expected to create changing pressure fields. These instabilities can break the osmotic pressure symmetry along the membrane, and after division the reactions defining the metabolism must be able to trigger a new growth-division cycle.

Let us first present the specific set of chemical reactions defining our basic metabolism. Several choices are possible and here we use one of the simplest scenarios found. As discussed by Morowitz (1992), the logic of replication could be separated from chemical constraints, but in order to be able to test the feasibility of a given minimal protocell model, we should consider chemically reasonable sets of reactions to be implemented. Here we present such a simple, but chemically reasonable, scenario. The set of reactions used here are(3.2)

(3.3)

(3.4)

In our model, reactions (3.1) and (3.2) can take place only inside the vesicle (for example due to the existence of inhibitory conditions outside it). Here *R*_{1} and *R*_{2} are the basic reagents, which are continuously and uniformly pumped from a source located at the limits of the system. The different substances involved in these reactions can cross the membrane with certain permeability and diffuse. The concentrations of the different molecules are denoted as *C*_{z}(*r*, *t*), where *z*∈{*R*_{1}, *R*_{2}, *g*_{1}, *g*_{2}, *g*_{0}}. For notational simplicity, the spatio-temporal dependence is not explicitly written. The following set of partial differential equations describes the dynamics of the proposed system:(3.5)(3.6)(3.7)(3.8)(3.9)Here, *ρ*_{1}, *ρ*_{2} and *ρ*_{3} are the degradation rates of *g*_{1}, *g*_{2} and *g*_{0}, respectively. The diffusion coefficients are denoted by , , , and . For simplicity, we assume that the value of the diffusion coefficients is the same inside and outside the membrane. Finally and are the constant rates of reagent supply.

This set of chemical reactions is able to trigger the emergence of a non-uniform spatial concentration of morphogens as a consequence of Turing-like instabilities. These instabilities are generated by the autocatalytic reactions (3.1) and (3.2) associated with the inhibitor effect of reaction (3.1). The previous model ((3.4) to (3.8)) can be numerically solved by using a spatial discretization of the surface domain *Γ*, and considering zero-flux boundary conditions at the limits *γ* of the domain, i.e. at the membrane.

We start with an initially homogeneous state where with a constant and *ξ*(*r*, 0) a small noise term (Meinhardt 1982). As a result of the previous set of interactions, concentrations change until they achieve a steady state. In figure 3 an example of the spatial distribution of *g*_{1} and *g*_{2} is shown, using the parameters given in table 1. As expected from a symmetry-breaking phenomenon, the two morphogens get distributed in separated (exclusive) spatial domains. Each one tends to concentrate in one of the poles. These effects, coupled with membrane growth, will be exploited to design an active mechanism for controlled membrane division.

A second component of our model involves membrane growth. The cell membrane will grow as a consequence of the continuous input of molecules or aggregates available from an external source. As a consequence of this process, the boundary *γ* (which now allows diffusion with the external environment) is not rigid anymore. At each time, *γ* will be a time-dependent function *γ*(*t*). The concentration at each instant depends on the number of molecules *n*_{j} and the volume *V*. This dependence can be indicated as follows:(3.10)Considering that *c*_{j}(*t*)=*n*_{j}(*t*)/*V*(*t*), equation (3.9) becomes(3.11)The first term on the right-hand side accounts for the change in concentrations associated with changes in the number of moles *n*_{j} inside the protocell, assuming constant volume. The second term accounts for the change in concentration due to changes in membrane volume associated with membrane growth. If the composition of the external solution does not change over time, the rate at which the externally provided compounds are incorporated into the membrane can be considered proportional to its area *A*:(3.12)where *T*_{d} is the time required until membrane size is duplicated.

Cell volume changes due to both net water flow and the growth of the membrane:(3.13)

In order to compute the net water in- and outflow, we must consider all the flows crossing the membrane. These flows are described by an additional set of equations. These equations account for the different interactions between all the elements of our system: water, solutes and membrane (Kedem & Katchalsky 1958; Patlak *et al*. 1963; Curry 1984). The net water flows can be expressed as follows:(3.14)where *L*_{p} is the hydraulic conductivity of the water; Δ*p* is the hydrostatic pressure difference between the interior and the exterior; *R* is the ideal gas constant; *T* is the temperature; and *σ*_{j} is the solute reflection coefficient for the *j*th substance (0 for a freely permeable solute, and 1 for a completely impermeable one).

Moreover, to properly compute the concentration changes through time in equation (3.10), we must also take into account the total substance flow inwards or outwards through the membrane. The flow can be expressed as (Curry 1984)(3.15)Here, *h*_{j} is the permeability of the *j*th substance (defined as the rate at which molecules cross the membrane), and(3.16)is the so-called Peclet number.

## 4. Membrane shape

Membrane growth, as described by equations (3.11) and (3.12) will be affected by the Turing instabilities generated by reactions (3.1) to (3.3). The membrane expansion is followed by a loss of spatial symmetry due to the effects of non-uniform osmotic pressures along the membrane surface *γ*(*t*).

The osmotic pressure at each membrane point is related to the current gradient of concentration between both membrane sides. At each point *r*∈*γ*(*t*), the osmotic pressure value generated by the *j*th substance can be written as(3.17)where *k*_{j} is a constant. For very low concentrations, we have *k*_{j}=*RT*, where *R* is the ideal gas constant and *T* is the temperature (if the concentrations are expressed in mols per liter). The osmotic pressure at one point *r* of the membrane and at time *t* can be calculated by adding the pressure generated by each substance as follows:(3.18)

Finally, we must take into account the contribution of the surface tension and the bending elasticity to the total pressure. This gives (for our two-dimensional system)(3.19)where *γ*_{0} is the surface tension coefficient, and *k* is the elastic bending coefficient. Here *R*_{s}(*r*) is the local radius of curvature, and *r*_{0} is the spontaneous radius of curvature.

For simplicity, we focus our attention on a two-dimensional model. The method employed to study the evolution of membrane shape has been previously presented in Macía & Solé (2006). This method considers only the local effect of osmotic pressures computed at each membrane point. Calculations are performed on a grid (Schaff *et al*. 2001). Figure 4 indicates how this discrete approximation is possible. In the lattice, there are three types of discrete elements: internal elements, which cover all the internal medium; external elements; and membrane elements, which cover the real membrane and are in contact with both internal and external elements.

Membrane shape can be described by a set of ‘characteristic points’ *Q*_{k}, each one associated with one membrane element. The position of each of these points can change dynamically as a consequence of pressure changes. In a first approximation, the displacement is proportional to the total pressure described by equation (3.18)(3.20)Here Δ*t* is the discrete time-interval used in the computation and *b* is a constant. This constant value cannot be arbitrary. The position of the characteristic points defines the membrane shape and size. Such shape and size must be in agreement with the membrane size and volume as determined by equations (3.11) and (3.12).

## 5. Results

The previous set of equations and boundary conditions allows the development of membrane growth and instability, using a realistic set of parameters (table 1). As discussed in §3, the system of chemical reactions generates steady non-uniform spatial concentrations. The question was to see if such spatial instabilities could trigger membrane changes leading to self-replication. In this section, we summarize our basic results showing that our model indeed displays the expected Turing-induced replication cycle.

As shown in figure 5, a spatial instability develops rapidly (figure 5*a*). The container where the metabolic reactions (3.1) to (3.3) are confined is a vesicle membrane. This membrane can be deformed and grow due to the incorporation of new, externally provided, precursor molecules. The growth in membrane size, given by equation (3.11), is related to the volume increase given by equation (3.12). When cell area increases, the internal volume also increases and the Turing instabilities move the maximal concentrations of *g*_{1} and *g*_{2} in opposite directions. Assuming that water and the different substances can cross the membrane with certain permeability, the non-uniform concentration distribution generates a non-uniform osmotic pressure along the membrane. Due to this pressure, the membrane can be deformed with the characteristic shape described in figure 1. In this context, the metabolic reactions (3.1) to (3.3) have an active control on the membrane growth and shape.

Figure 5 shows the evolution of the concentration profiles of *g*_{1} and *g*_{2} coupled with the membrane expansion process. As cell volume increases (in our two-dimensional model this is represented by the internal area |*Γ*|), it enhances the spatial segregation of morphogens, due to the increase in the size of the domain. In fact, in these regions where the concentrations of *g*_{1} and *g*_{2} are maximal there is a maximal expansive osmotic pressure. This pressure enhances the expansion in the poles and, as a consequence, the compression in the equatorial zone. The expansion taking place in the two poles enhances the separation of *g*_{1} and *g*_{2}, and so on. The coupling of these two mechanisms, each one enhancing the other, creates a controlled membrane deformation. Beyond a certain critical point membrane splitting becomes energetically favourable. In figure 6, we show the pressure distribution along the membrane in different simulation steps, consistently matching the theoretical pattern indicated in figure 1. Eventually, the narrow membrane division is a singularity and must be specifically introduced in the simulations (Macía & Solé 2006).

After division, concentration distributions in each daughter cell are not the same as those at the start of division cycle. Figure 5*a* shows the initial steps of the first cycle division, where the maximal concentration of *g*_{1} and *g*_{2} are comparable. However, figure 5*d* shows the situation right after division. In each cell there is a clearly dominant substance. Figure 7 shows the evolution on one daughter cell (the left one of figure 5*d*). As shown in figure 7*a*–*c* the minority substance regenerates and a new division cycle takes place. This is possible due the fact that, in spite of the kinetically symmetric features of both *g*_{1} and *g*_{2} (table 1), their respective reagents are different from the reagent of reaction (3.2).

The concentration of the dominant morphogen (*g*_{2} in figure 7*a*) has a growth limitation due to the substrate (*R*_{2}) depletion. This depletion is a consequence of the flow penetration limitation imposed by the membrane permeability. However, for the minority substance (*g*_{1}), the substrate consumption (*R*_{1}) is lower due to the low concentration of the autocatalytic substance *g*_{1}. At the beginning, the inflow rate of *R*_{1} is greater than its consumption. If the inhibitory effect of reaction (3.3) is not enough to limit the autocatalytic growth of *g*_{1}, the concentration can grow until the consumption of substrate *R*_{1} is faster than the inflow rate. Then the depletion of substrate *R*_{1} becomes the dominant mechanism and the growth of *g*_{1} is limited. To accomplish these effects, it is required that the diffusion coefficients of *R*_{1} and *R*_{2} are larger than the diffusion coefficient of *g*_{1} and *g*_{2}.

These are common characteristics of Turing pattern formation in a wide range of scenarios (Meinhardt 1982). Furthermore, the permeability of *R*_{1} and *R*_{2} must be bigger than the permeability of *g*_{1} and *g*_{2} in order to accomplish a significant effect of osmotic pressure in the poles. Finally, the kinetic constants *k*_{1} and *k*_{2} must be greater than *k*_{3} to ensure that the inhibitory effect of reaction (3.3) enhances the substances separation without preventing the growth of the minority substance (table 1).

Our model assumes that the membrane remains continuously closed through time, as it grows following equation (3.11).1 In the model, the position of the characteristic points *Q*_{k} determines the size and shape of the membrane. In order to have a physically consistent simulation, the size of the membrane calculated with equation (3.11) must be in agreement with the size as determined by the spatial location of the characteristic points. Figure 7*d*,*e* shows the dynamics of cell growth and replication in terms of membrane size. In figure 7*e*, we can see the agreement between membrane size as calculated by equation (3.11) and the one derived from the characteristic point locations, as determined by equation (3.19). In figure 7*d*, we display the membrane size evolution along three division cycles. The small differences between the consecutive cycles are an artefact of the model discretization.

## 6. Discussion

We have presented a mesoscopic, minimal cellular model defined in terms of a closed membrane with a simple internal metabolism. These metabolic reactions are able to create Turing-like instabilities. Specifically, a mechanism leading to lateral inhibition associated with exclusive states has been used (Meinhardt 1982). These instabilities, coupled with membrane growth, provide an active method for controlled membrane deformation and have been shown to trigger cell replication. The basic mechanism is related to a non-uniform osmotic pressure distribution along the membrane. Although a specific mechanism has been presented here, we have found other possible (more sophisticated) scenarios where this also occurs (J. Macia & R. V. Solé 2006, unpublished data).

Since spatial instabilities play an important role in many natural processes, the design of mechanisms based on these instabilities could be relevant to the synthesis of artificial protocells and even for understanding prebiotic scenarios of cellular evolution, where the sophisticated division mechanisms of the current cells were not present. In this context, the set of metabolic reactions can be arbitrarily generalized, thus opening the door for many different types of membrane–metabolism couplings. Future work should explore the possible paths towards the experimental synthesis of these protocells from available molecular structures, their potential evolvability and further extensions to more complex metabolic networks (Kaneko & Yomo 2002).

## Acknowledgments

The authors thank the members of the Complex System Lab for useful discussions. Special thanks to Carlos Rodriguez-Caso for helpful comments and suggestions. This work has been supported by EU PACE grant within the 6th Framework Program under contract FP6-0022035 (Programmable Artificial Cell Evolution), by McyT grant FIS2004-05422 and by the Santa Fe Institute.

## Footnotes

↵1 Other parameters can be unable to keep the membrane growing or, instead, make it grow too fast. In those cases, membrane breaking can occur.

- © 2007 The Royal Society