Stochastic simulations of minimal cells: the Ribocell model
 Fabio Mavelli^{1}Email author
Affiliated with
DOI: 10.1186/1471210513S4S10
© Mavelli licensee BioMed Central Ltd. 2012
Published: 28 March 2012
Advertisement
DOI: 10.1186/1471210513S4S10
© Mavelli licensee BioMed Central Ltd. 2012
Published: 28 March 2012
Over the last two decades, lipid compartments (liposomes, lipidcoated droplets) have been extensively used as in vitro "minimal" cell models. In particular, simple and complex biomolecular reactions have been carried out inside these selfassembled micro and nanosized compartments, leading to the synthesis of RNA and functional proteins inside liposomes. Despite this experimental progress, a detailed physical understanding of the underlying dynamics is missing. In particular, the combination of solute compartmentalization, reactivity and stochastic effects has not yet been clarified. A combination of experimental and computational approaches can reveal interesting mechanisms governing the behavior of micro compartmentalized systems, in particular by highlighting the intrinsic stochastic diversity within a population of "synthetic cells".
In this context, we have developed a computational platform called ENVIRONMENT suitable for studying the stochastic time evolution of reacting lipid compartments. This software  which implements a Gillespie Algorithm  is an improvement over a previous program that simulated the stochastic time evolution of homogeneous, fixedvolume, chemically reacting systems, extending it to more general conditions in which a collection of similar such systems interact and change over the course of time. In particular, our approach is focused on elucidating the role of randomness in the time behavior of chemically reacting lipid compartments, such as micelles, vesicles or micro emulsions, in regimes where random fluctuations due to the stochastic nature of reacting events can lead an open system towards unexpected time evolutions.
This paper analyses the socalled Ribocell (RNAbased cell) model. It consists in a hypothetical minimal cell based on a selfreplicating minimum RNA genome coupled with a selfreproducing lipid vesicle compartment. This model assumes the existence of two ribozymes, one able to catalyze the conversion of molecular precursors into lipids and the second able to replicate RNA strands. The aim of this contribution is to explore the feasibility of this hypothetical minimal cell. By deterministic kinetic analysis, the best external conditions to observe synchronization between genome selfreplication and vesicle membrane reproduction are determined, while its robustness to random fluctuations is investigated using stochastic simulations, and then discussed.
In recent years, many researchers have been actively working in the field of the de novo synthesis of the artificial cell [1–5], i.e. a cell made from scratch using both synthetic and natural compounds. This scientific challenge has many relevant aspects: first of all, it can reinforce the theory of abiogenesis in the origins of life debate [6], proving that life can emerge spontaneously in a test tube, at least in suitable experimental conditions. Furthermore, the possibility to create a population of artificial cells programmed for the synthesis of chemical compounds of pharmacological and industrial interest is a significant biotechnological goal [7, 8]. Artificial cells can also be envisaged, in the maybe not too distant future, as microscopic diagnostic and pharmacological labs to be delivered into the human body in order to synthetize and release drugs as a response to an external stimulus in presence of a disease [9].
The topic of artificial cells is strongly related to the minimal cell notion that defines the simplest cell to be considered alive and then to be experimentally implemented. A minimal living cell, or protocell, can be defined as the minimum supramolecular bounded structure based on the lowest number of molecular species and metabolic processes that is capable of selfmaintaining, selfreproducing and evolving [10]. SelfMaintenance is the necessary condition that must be fulfilled by the protocell to be considered alive, i.e. it must stay in a steady state where all its constituents are continuously synthesized and refurbished [11]. Nevertheless, this is not sufficient to implement cellular life as we actually know it; in fact, real cells are also able to grow and to selfreproduce. On the other hand, SelfReproduction can be seen as a consequence of the cellular metabolism that can keep the minimal cell in a stationary state or bring it to a continuous growing regime that forces the organism to divide in order to maintain its internal coherence, i.e. its stability. Both these features are individual properties of a single cell that can be observed during its life time. In contrast, Evolvability is a collective property that can be exhibited only by a population of cells and on a time scale of several generations, according to a Darwinian selection mechanism [6, 12].
Some years ago, Szostak and colleagues proposed a minimal cell prototype called the Ribocell: the RNAbased cell, that in principle can exhibit all the three properties to be considered a living protocell [13]. This theoretical cellular model consists in a selfreplicating minimum genome coupled with the selfreproduction of the lipid vesicular container [14]. These authors envisaged the existence of two hypothetical ribozymes [13], one (R_{L}) able to catalyze the conversion of molecular precursors (P) into membrane lipids (L) and the other one (R_{P}) able to duplicate RNA strands. Therefore, in an environment rich in both lipid precursors (P) and activated nucleotides (NTPs), the Ribocell can selfreproduce if both processes, i.e. genome selfreplication and membrane reproduction (growth and division), are somehow synchronized.
Kinetic Parameters for the in silico Ribocell model at room temperature.
Kinetic Parameters  Values  Process Description  Ref. 

k_{SS}[s^{1}M^{1}]  8.8·10^{6}  Formation of dimers R_{c}R_{P} and R_{c}R_{L}  [28] 
k_{S}[s^{1}]  2.2·10^{6}  Dissociation of dimers R_{c}R_{P} and R_{c}R_{L}  [28] 
k_{R@S}[s^{1}M^{1}]  5.32·10^{5}  Formation of R@S  [29] 
k_{R@SS}[s^{1}]  9.9·10^{3}  Dissociation of Complexes R@S_{c}S  [29] 
k_{NTP}[s^{1}M^{1}]  0.113  Nucleotide Polymerization in Oleic Vesicle  [31] 
k_{L} [s^{1}M^{1}]  1.7·10^{3}  Lipid Precursor Conversion*  [30] 
k_{in} [dm^{2}s^{1}]  7.6·10^{19}  Oleic acid association to the membrane  [26] 
k_{out} [dm^{2}s^{1}]  7.6·10^{2}  Oleic acid release from the membrane  [26] 
P_{P} [cm·s^{1}]  4.2 10^{9}  Membrane Permeability to Lipid Precursor  
P_{NTP} [cm·s^{1}]  1.9 10^{11}  Membrane Permeability to Nucleotides  [31] 
P_{W} = P_{S}  0.0  Membrane Permeability to W and genetic staff  
P_{aq}[cm·s^{1}]  1.0·10^{3}  Oleic Acid Membrane Permeability to Water  [32] 
In this paper, we first apply this model to 100baselong ribozymes, in an attempt to find the best experimental conditions to reduce so as Ribocell life time. By using the deterministic approach, the robustness of the stationary growth and division regime will be investigated in terms of external substrate concentrations, vesicle size and initial ribozyme amount in order to define optimal external conditions for Ribocell selfreproduction.
Therefore, the influence of ribozyme length will also be explored in the optimal external conditions by ranging strand size from 20 to 200 bases in length and keeping all the other kinetic parameters constant. 20 bases is in fact the minimum length required to observe a folded RNA structures, i.e. a structure that can reasonably exhibit catalytic action. On the other hand, entities of about 200 nucleotides have been suggested as plausible ancient protoribosomes [18] even though, more recently, smaller subunits of 60 nucleotides have also been considered as plausible candidates [19]. Moreover, the dependence of Ribocell time behavior on the kinetic constants of RNA dimer formation and dissociation will also be studied.
Finally, stochastic simulations will be performed in order to test the robustness of the ribocell base on 100base length ribozymes in optimal external conditions, with the aim of elucidating the role of intrinsic and extrinsic stochasticity on the time behavior of a protocell population.
Selfreproducing vesicles are compartmentalized chemicallyreacting systems where selfassembly processes are coupled with chemical reactions that produce amphiphilic molecules. To study the time behavior of these systems, we use both a deterministic and a stochastic approach in order to get insights into the average behavior of the protocell population and, at the same time, to elucidate the role of random fluctuations. Given a certain minimal cell model, i.e. a reaction mechanism with all the required parameters (kinetic constants, permeability coefficients, initial concentrations), the deterministic analysis can be done by numerically solving the ordinary differential equation set (ODES) or by analytically integrating an approximated reduced set of differential equations. Examples of the latter approach can be found in our previous works where selfreproducing micelles [20] and vesicles [21] were studied. In order to take into account the stochastic effects, some years ago we developed a Monte Carlo program [22] based on the Gillespie's direct method [23] designed for coupling amphiphile selfassembly [24] with chemical reactions and hydrophobic solute absorption [25] in a homogeneous well stirred macroscopic reactor. More recently, this program has been improved and a computational platform called ENVIRONMENT [26] was later developed to cope with the more general case of a collection of reacting lipid compartments that interact and change over time. In the next subsections, the general model for a selfreproducing vesicle will be recalled and discussed [17].
where n _{i} ^{C} are the molecular numbers of species X_{i} (i = 1,2 ... N) present in the vesicle aqueous core and n _{L} ^{μ} is the number of amphiphiles X_{ L } (1≤ L ≤N) in the membrane. V _{C} is the water internal volume. In the stochastic approach, all n _{i} ^{C} and n _{L} ^{μ} are discrete integer numbers and there exist as many arrays as vesicles in the systems, while in the deterministic analysis there is a single array with real values that represents the average time state of the entire reacting vesicle population.
Deterministic rates against propensity density probability for reacting and transport events.
Event  Deterministic Rate (Ms)^{1}  Propensity Density Probability s ^{1} 

Internal Chemical reactions ^{(a)}



Solute X_{n} membrane transport ^{(b)} 


Membrane Lipid Release 


Membrane Lipid Uptake 


Water Flux ^{(d)} 


ε and η being the osmotic and dividing tolerance coefficients, respectively. In fact, vesicles in hypotonic solutions can swell, stretching the membrane until they reach a critical state: ϕ = 1ε. The osmotic tolerance ε can be experimentally determined by measuring the maximum difference in osmolite aqueous concentrations between the internal core (C_{ T }^{ C }) and the external environment (C_{ T }^{ E }) that vesicles can bear. For oleic acid, ε was found equal to 0.21 [27]. In our model, at the critical bursting point, vesicles are assumed simply to break down, releasing all their internal content into the external environment and remaining in solution as flat bilayers, since bilayer sealing processes are not considered in this model. On the other hand, deflated vesicles are supposed to be able to divide in order to minimize the bending energy. The dividing condition is reached when they can form two equal volume spherical daughters (ϕ = 2^{1/3}). So η introduces a tolerance that is linked to the relative flexibility exhibited by any membrane. However, as a simplifying assumption in this work, η = 0 will be supposed. As already mentioned, after each division, all the molecular content of the mother vesicle is halved in the deterministic approach, while it is randomly distributed between the twin daughters in the stochastic simulations.
Figure 1 reports the internal metabolism proposed for the Ribocell. Both pairs of RNA strands reversibly associate (1) and these processes are shifted towards the dimer formation and are strongly dependent on temperature. The replication of any RNA strand is catalyzed by the polymerase R_{P} according to the steps in bracket (2). This process is described as a catalytic templatedirected addition of mononucleotides with high fidelity and processivity. It starts with R_{P} binding any of the monomeric strands S (S = R_{P}, _{c}R_{P}, R_{L} and _{c}R_{L}) to form the complex R@S. This complex will then initiate the polymerization of the conjugate strand _{c}S, by coupling and iteratively binding the complementary bases and releasing the byproduct W. When the strand _{c}S has been completely formed, the polymerase ribozyme releases the new dimer. Finally, the ribozyme R_{L} catalyzes the conversion of the precursor P into the lipid L (3).
Table 1 shows the values of all kinetic constants and membrane permeabilities used in this work, with their respective references. In this work, the length of the two ribozymes is assumed to be 100 nucleotides long (20 being the minimum base number for observing a folded RNA conformation). Both R_{L} and R_{P} are created with a random sequence of bases, and they are assumed to have similar kinetic behaviors for the sake of simplicity. The kinetic constants of formation k _{SS} and dissociation k _{S} of both dimers were set equal to the values experimentally measured for a sequence of 10 nucleotides [28]. The kinetic constants for both complex formation R@S and complex dissociation R@S_{c}S (S = R_{P}, _{c}R_{P}, R_{L} and _{c}R_{L}) were set equal to those measured for the human enzyme βpolymerase [29] that binds a DNA single strand and dissociates from a DNA double helix, respectively. The rate constant for the catalytic synthesis of lipid k _{ L } was taken to be 10^{5} times larger than that of the splicing reaction, catalyzed by the hammerhead ribozyme [30]. The kinetic behavior of different nucleotides is assumed to be the same and a single value is assigned to k _{NTP} derived from experimental data simulations (De Frenza private communication) of the DNA template directed synthesis in fatty acid mixed vesicles [31]. In the same way, the common value of the membrane permeability to activated nucleotides was also estimated, while the membrane is assumed impermeable to genetic material. The kinetic constants of the membrane/aqueous solution lipid exchange for oleic acid vesicles are taken from a previous work [26] where they were obtained by simulating the competition between isotonic and osmotically swollen oleic vesicles [27]. The amphiphile head area α _{ L } = 0.3 nm^{2} and the osmotic tolerance ε = 0.21 are defined according to data reported in literature for oleic acid vesicles [27]. The only two parameters assigned arbitrarily are therefore the membrane permeability to the byproduct: P _{W} = 0.0 cm/s, based on the assumption that W is a charged species, and the permeability to the precursor: P _{P} = 0.42·10^{8} cm/s, corresponding to the oleic acid membrane's permeability to Arabitol [32] and comparable to those of similar organic compounds.
A common simplifying assumption to both approaches is to consider the external concentrations of nucleotides (NTPs), lipid precursor (P) and inert compound (B) to be constant throughout the timecourseof the process, thanks to an incoming flux of material in the reactor vessel: continuous stirred tank reactor approximation. The byproduct (W) concentrations is also assumed to be constantly equal to zero outside. For all these compounds, except B, the internal aqueous concentration is zero at the beginning. Moreover, all the calculations are performed starting from an initial isotonic condition, so that the internal concentration of B is properly adjusted in order to counterbalance the presence of the lipid precursor and nucleotides outside and the genetic material inside, respectively.
As has been pointed out before, given the set of kinetic parameters reported in Table 1 we are looking for a set of initial concentrations that can allow the Ribocell to reach a stationary regime of growth and division, i.e. a dynamic state where the protocell grows and divides producing two twin daughters the same size as their mother at the beginning of its life cycle.
it is easy to show(unpublished observations) for growing protocells that synchronization can take place only if γ = 1, while if γ > 1 the volume will increase much faster than the membrane surface and the vesicle can become energetically unstable, leading to an osmotic burst. On the other hand, if γ < 1 then the vesicle will divide decreasing in size generation by generation. It is beyond the scope of this work to go into the mathematical details of this formula, so the interested reader should refer to an incoming paperfor a more detailed discussion.
Nevertheless, it is worthwhile to keep in mind that v _{S} is essentially proportional to the rate of lipid synthesis, since amphiphile uptake by the membrane is very fast when the concentration of lipids is above the equilibrium value. Instead, having assumed the external value C _{T} ^{E} to be constant and a vesicle being in a osmotic balanced condition: C _{T} ^{E}≈ C _{T} ^{C} = V _{C} N _{T} ^{C},v _{V} is driven by the overall internal N _{T} ^{C} population rise. Therefore, since N _{T} ^{C} increases essentially owing to the waste production that takes place with ribozyme selfreplication and lipid synthesis, this has the effect of coupling the membrane reproduction with the genome replication and with the volume growth: osmotic synchronization.
In the present paper, we firstly explore the dependence of the stationary regime on the external concentrations of substrates for ribozymes 100 nucleotides long setting k _{L} to 1.7·10^{3}s^{1}M^{1}, i.e. the minimum value in the previously observed synchronization range. The aim of this preliminary deterministic study is to find the optimal initial conditions in order to achieve a stationary regime with the shortest life time. All the outcomes are reported in additional file 1.
Figure 3. Two equal fractions of R_{L} and _{c}R_{L} strands are present, both around 33%, while the percentages for R_{P}, and _{c}R_{P} are lower: 24% and 10%, respectively. Due to the high k _{ss} value, the ribozymes are mainly present inthe form of dimers and this accounts for the equal fractions observed for the lipase ribozymes, while the percentage of R_{P} is greater than that of _{c}R_{P}, since some R_{P} strands are involved as catalyzers in template duplication. This also explains why the overall percentage of polymerase ribozymes (~34%) is lower than that of lipase ribozymes (~66%) in fact, not all polymerase ribozymes are available as templates for duplication. Moreover, this asymmetry is amplified as long as the total concentration of the genetic material increases, see data in additional file 1.
Stochastic simulations were performed by means of the parallel version of ENVIRONMENT, running 32 statistically equivalent simulations of a 10ribocell solution on different CPUs. Therefore, the outcomes were obtained as averages from a population of 320 vesicles. Kinetic parameters used for simulations are those reported in Table 1 while the initial conditions are shown in bold by additional file 1. At each cell division, only one of the two offspring was kept while the other was discarded in order to reduce computation time, thus keeping the number of monitored vesicles constant. This is in agreement with the assumption that the external concentrations of all substrates are fixed due to an incoming flux of material, i.e. the substrates cannot ever be exhausted.
In this paper, we applied an already published Ribocell in silico model [15, 16] to the case of two hypothetical ribozymes 100 nucleotides long by using the kinetic parameters reported in Table 2. The kinetic constant for the lipid formation, k _{L}, was set equal to 1.7E+3s^{1}M^{1}, the lowest value that exhibited a stationary regime of growth and division in a previous work [16] in which ribozymes 20 nucleotide long were assumed. Length of 100 nucleotides was chosen as a compromise between the need to reduce calculation times and the choice of a plausible ribozyme size, keeping in mind that the RNA subunits in ribosomes devoted to protein elongation have a comparable length [18, 19].
By means of deterministic analysis, the robustness of the stationary regime was also investigated as a function of the initial conditions, the length of ribozymes and the kinetic constants of the RNA dimerization. For 100base long ribozymes, the best experimental conditions in terms of the external concentrations have been found: [N _{ex}] = [P _{ex}] = 1.0E2M and [I _{ex}] = 0.3M in order to observe a stationary regime with the lowest division time: Δt _{20} = 68.3 days. A so high protocell life time can be mainly ascribed to the RNA reversible association that is shifted towards the dimer formation rendering the concentration of the catalyzers R_{L} and R_{P} very low. This has been confirmed by analyzing the dependence of the life cycle on the thermodynamic constant of ribozyme dimerization. A small influence on the stationary regime was observed on changing the length of the RNA strands: an increase in filament size determines higher life times and a lower amount of genetic material. This effect is more pronounced when the R_{P} length changes. On the other hand, the external concentrations of substrates and inert compounds appear to highly affect the stationary regime in terms of vesicle size, genetic material amount and genome composition. Moreover, deterministic calculations have also shown that the stationary regime can be reached from very different initial genome composition, although when too much free R_{L} ribozymes are present at the beginning the death for ribozymes segregation can be observed since the protocell divides too quickly before the genome replication. Stochastic simulations have been done starting from a population of 320 identical ribocells with an initial genome composed by N _{0} = 1,10 and 100 dimers of both R_{c}R_{L} and R_{c}R_{P}. The analysis of simulations outcomes shows that the ribocell time behaviors is highly influenced by random fluctuations. Since the genetic material is randomly distributed at each cell division, this can produce different type of protocells, ranging from empty vesicles to genuine ribocells, their internal metabolism being highly influenced by the presence of the catalytic RNA strands. In fact, deterministic analysis cannot take into account the disappearance of ribozymes due to a vesicle division, since this approach simply halves the genetic amount and follows the reacting molecule time courses in terms of population averages, i.e. real positive numbers that can be less than one without being zero. On the other hand, the stochastic simulations are more realistic to the random loss of ribozymes from the genome being capable of describing a population of protocells with completely different time behaviors. As a consequence, the simulation outcomes show that ribocells are not enough robust to survive to random fluctuations. In fact only about the 5% of the initial population survive as genuine ribocells after 1525 generations and on a longer time window they are destined for extinction. Furthermore, the time course of each single protocell is also greatly influenced by intrinsic stochasticity in particular by the time fluctuations of the RNA dimer dissociation. In fact, when all the RNA strands are associated in dimers, protocells remain in a lazy phase, whereas free R_{L} monomers induce fast growth and division steps and free R_{P} cause the fast RNA replication without changing the vesicle size appreciably. Therefore these two process are synchronized only by chance and this also represents a reason of weakness of this model protocell.
In order to implement experimentally ribozymesbased minimal cells two main improvements are necessary. As first, more free monomers of both R_{L} and R_{P} must be available in the vesicle core so that the ribocell life cycle will be speeded up and the division time lowered. This can be achieved by increasing the working temperature since it has been recently show that fatty acid vesicles are stable up to 90°C [33] and the efficiency of the selfcatalyzed replication of RNA strands increases with a temperature rise. This is also in agreement with our results since the k _{SS} value used can be considered as the appropriate 100base long RNA association constant for a higher temperature than 25°C [34]. For a more detailed theoretical analysis at high temperature, it is necessary, of course, to estimate the kinetic constants for all the steps involved in the internal metabolism. This will be the topic of a future work. The second necessary improvement consists in finding a way to really synchronize the genome selfreplication and the membrane reproduction. This is a much more complex task to achieve. Working with high concentration of the genetic material, it can avoid, or at least reduce, the ribozymes segregation and this should be compatible with a high working temperature. The best strategy could be to have a fine control of the R_{c}R_{L} dissociation since when free R_{L} monomers are present in the aqueous core the membrane growth quickly and the division takes place very soon. Thus the R_{c}R_{L} dissociation can act as a trigger for the membrane growth and division.
Finally rephrasing the George Box famous sentence, we are aware that this in silico Ribocell model is in a some way wrong, but we hope it might inspire researchers involved in the lab implementation of the ribozymesbased minimal cell.
This work has been funded by the MIUR (PRIN2008FY7RJ4_002). Parallel calculations were done at CINECA Super Computing Italian Centre thanks to the approved ISCRAC project (HP10CVJLGZ). The authors would like to thank Anthony Green for his careful review of the English in this manuscript.
This article has been published as part of BMC Bioinformatics Volume 13 Supplement 4, 2012: Italian Bioinformatics Society (BITS): Annual Meeting 2011. The full contents of the supplement are available online at http://www.biomedcentral.com/14712105/13/S4.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.