Abstract
Recent experimental work in the field of synthetic protocell biology has shown that prebiotic vesicles are able to ‘steal’ lipids from each other. This phenomenon is driven purely by asymmetries in the physical state or composition of the vesicle membranes, and, when lipid resource is limited, translates directly into competition amongst the vesicles. Such a scenario is interesting from an origins of life perspective because a rudimentary form of celllevel selection emerges. To sharpen intuition about possible mechanisms underlying this behaviour, experimental work must be complemented with theoretical modelling. The aim of this paper is to provide a coarsegrain mathematical model of protocell lipid competition. Our model is capable of reproducing, often quantitatively, results from core experimental papers that reported distinct types vesicle competition. Additionally, we make some predictions untested in the lab and develop a general numerical method for quickly solving the equilibrium point of a model vesicle population.
Introduction
A fundamental problem in biology concerns the origins of an innovation that allowed the development of organisms in our biosphere, beyond complex chemical reaction networks: the emergence of cells^{1,2}. Cells define a clear scale of organization and, given their spatially confined structure, they constitute efficient units where molecules can easily interact, coordinate their dynamical patterns and establish a new level of selection. Although it is often assumed that there was a transition from some type of ‘lessorganised’ prebiotic chemistry (probably including catalytic cycles) to a cellbased living chemistry, little is yet known about the potential pathways that could be followed to cross it. Once in place, protocell assemblies would require available resources for their maintenance and, thus, would naturally get inserted in diverse competitive dynamics in which the main selective unit would be the whole protocellular system. In this context, aggregatelevel evolution is the right scale of analysis to be considered.
Different types of protocellular systems of diverse complexity have been studied from a theoretical standpoint^{3,4,5,6,7,8,9,10,11}. In particular, by considering the coupling of a template carrying information with vesicle replication and metabolism, it has been shown that Darwinian selection is the expected outcome of competition in a protocellular world^{12}. In a more simple scenario for autopoietic (i.e. selfproducing) vesicles in a homeostatic regime, previous numerical simulations suggest that random fluctuations can also act as ‘selection rules’ for the more robust individuals^{13}. Early preDarwinian stages in the development of biological organisms in which supramolecular systems could still be disconnected from information (i.e. closer to elementary forms of metabolism and strongly constrained by the molecular diversity of the available chemical repertoire) ought to be further explored. What type of competition and cooperation processes were at work in the chemical world leading to the emergence of early protocells? Processes able to favour asymmetries in the chemical composition of vesicles should be expected to play a relevant role in this context, defining the conditions under which protocellular assemblies could thrive.
Recent laboratory experiments have actually demonstrated how differences in the composition or physical state of the vesicle membrane can drive competition for simple amphiphilic molecules (typically fatty acids), present in solution as free monomers at low concentration. First, Cheng and Luisi^{14} observed competition between pure oleic acid and POPC vesicles, where each of these vesicle populations had different initial size distributions. In all the studied cases, the final size distribution was found near to the initial one of the POPC vesicles, suggesting that oleic acid molecules were rapidly absorbed by POPC aggregates. Then, Chen et al.^{15} reported competitive dynamics in a population of fatty acid vesicles, whereby vesicles that were osmotically swollen by an encapsulated cargo of RNA (or sucrose) stole lipids from their empty, osmotically relaxed counterparts by virtue of absorbing monomers more quickly from the solution. They studied both oleic acid and POPC vesicles, but only in the former case was competition observed. This distinctive behaviour of fatty acid vesicles has been theoretically rationalised^{16} by assuming that double chain phospholipids are taken up from solution by the vesicle membranes five orders of magnitude more slowly than single chain fatty acid molecules.
More recent experimental work has turned attention to other possible selective advantages of protocells, such as phospholipid^{17} and peptide^{18} driven competition amongst vesicles. Instead of membrane tension, the main factor for competition here is the different composition of the membrane; singlechain fatty acids are mixed with double chain amphiphiles or with a different type of surfactant molecule, like sufficiently hydrophobic peptides. In the case of phospholipiddriven competition, oleic acid vesicles endowed with a membrane fraction of phospholipid are observed to take fatty acid molecules from phospholipiddeficient neighbours, who shrink, whilst the former grow and keep their potential for reproduction.
In this paper we develop a mathematical model of a competing population of vesicles, with the motivation to explore and test possible mechanisms underlying lipid competition phenomena. The model is based at the coarsegrain level of lipid kinetics, following the approach of Mavelli and RuizMirazo^{16}. Using physically realistic parameters such as lipid molecule sizes, vesicle aggregation numbers and critical vesicle concentrations (CVC) as detailed in Table 1, we are able to qualitatively and often quantitatively reproduce results from two key experimental papers describing phospholipiddriven^{17} and osmoticallydriven^{15} competition.
In the model, a vesicle in a population absorbs and releases amphiphiles to and from its membrane at rates that depend on the current physical properties of that particular vesicle (such as membrane composition or extent of osmotic tension). To take account of phospholipiddriven competition, we build into the lipid kinetics two basic physical mechanisms, which have been postulated in the literature to underlie asymmetric growth dynamics in this context: the indirect effect and the direct effect, as will be named in this work. The first one refers to the decrease of amphiphile release processes simply due to the fact that other surfactant molecules are present in the membrane and the second to the immediate influence that these surfactant molecules could have on the amphiphiles (see Fig. 1).
More generally, this work forms part of our endeavour to try to develop a formalism that grasps the lipid kinetics involved in vesicle selfassembly under controlled conditions (pH, temperature, etc.). In contrast to the kinetics of chemical reaction networks which have been extensively modelled by the Mass Action Kinetics (MAK) and Stochastic frameworks^{19,20}, membrane lipid kinetics have been largely underexplored in the literature, due to the inherent complexity of supramolecular structures. Nevertheless, models coupling together membrane and metabolism kinetics will be a crucial cornerstone in order to build a systems understanding of the dynamic properties and organization of protocells, ultimately biological cells as well.
The paper is organised as follows. The remainder of the introduction serves to both introduce our kinetic model in detail and to perform a meanfield analysis of it. This analysis gives insight into why we should expect phospholipiddriven competition to result in a simplified version of our model. Then, the Results section summarises how well numeric simulations of the full kinetic model are able to reproduce experimental results and observations, including also some predictions for still untested protocell competition scenarios. In the Discussion section, we comment on some assumptions and other aspects of our approach and conclude the study. The Methods section at the end of the paper describes a fast numerical method for solving the final equilibrium state of the full model. This method was essential in producing the results figures in the paper. The Supplementary Material (online) explains aspects in more detail, including justification for some modelling choices and the vesicle mixing procedure assumed in order to compare our model with experimental observations.
Theoretical model of vesicle competition
The competition model (Fig. 2) involves a set of n vesicles , each one characterized by a quintuple and all embedded in a finite volume environment defined by a triple (Ω_{e}, L_{e}, B_{e}).
Each competing vesicle consists of a unilamellar (single bilayer) membrane of up to two different lipid types: singlechain fatty acid lipids (e.g. oleic acid, OA) and possibly a fixed number of doublechain phospholipids (e.g. dioleoylphosphatidic acid, DOPA). Membrane thickness is considered negligible and surface area of a vesicle, referred to as , is the waterexposed area of either of the bilayer leaflets. The L type lipids in the bilayer continuously exchange with the vesicle internal water pool and , whereas the phospholipids P are considered approximately stationary in the bilayer due to their comparatively slow exchange rate, in agreement with previous reported work using POPC vesicles^{16}. The internal water pool of each vesicle is considered a wellmixed chemical domain of volume Ω_{j} and hosts lipid monomers and also buffer species. Buffer species cannot permeate the bilayer but provide osmotic stability and they are also present in with constant number B_{e}.
Vesicles compete with each other by consequence of uptaking/releasing fatty acid monomers L from/to , which is a common limited resource. The initial system of vesicles is taken to be the result of mixing different vesicle populations and is a closed system in a nonequilibrium state. The system equilibrates to a final state following the dynamics described below, with some vesicles growing bigger in surface at the expense of others, which shrink. We ignore spatial correlations and the possibility of direct vesiclevesicle interactions and assume a wellmixed set of vesicles.
More precisely, each vesicle is considered to release lipids to both aqueous phases (at each side of the bilayer) at the equal rate of and absorb lipids from each phase at rate , where [L] is the molar concentration of lipid monomer in the respective phase. Functions r and u are defined later.
The uptake and release kinetics are symmetric on each side of the bilayer, which means that the lipid monomer concentration inside and outside each vesicle will be equal at equilibrium. Flipflop of the fatty acid L between membrane leaflets is considered very fast with respect to its uptake and release rates and thus a bilayer is modelled as a single oily phase; this simplification is supported by experimental work from Hamilton's lab^{21,22}. Conceivably, leaflet asymmetries could be created by the fact that the flipflop of deprotonated and protonated fatty acid molecules is not the same^{23}. However, such effects are considered of secondary importance and are disregarded in the present work.
Explaining the choice of L release kinetics, each fatty acid in a pure L membrane is considered to have a uniform probability per unit time k_{out} of disassociating from the membrane^{16}, while function r has been introduced in this work to take into account the direct effect. This function (0 ≤ r(ρ) ≤ 1) modifies the fatty acid release probability, based on the current molecular fraction of phospholipid P in a membrane . It is monotonically decreasing with increasing ρ, meaning that increasing phospholipid fraction generally decreases bilayer fluidity, slowing down the rate of L release from the membrane^{17}. In a first approximation, r was assumed linear:
where parameter 0 ≤ d ≤ 1 tunes how the lipid release rate is affected by phospholipid content (1 being maximally affected and 0 being not at all).
Conversely, lipid uptake kinetics reflect that the probability of uptaking a lipid L to the membrane is proportional to the density of lipid monomer in the immediate vicinity of the respective bilayer surface (i.e. the concentration of lipid in the surrounding medium), the area of surface available for absorption S_{μ} and function u, based on the dimensionless reduced surface The reduced surface encodes the surface area to volume ratio of a vesicle, when the latter internal volume is considered as a sphere: Φ = 1 denotes a vesicle perfectly spherical in shape, whereas Φ < 1 or Φ > 1 indicates a vesicle in osmotic tension or deflated respectively. Taking this into consideration, we define u as the following conditional function
to denote that lipid uptake is only increased when the the bilayer is stressed^{24}. Flaccid vesicles do not have extra enhancement of lipid uptake rate. Rationale for this function originated in the theoretical modelling of osmoticallydriven competition dynamics between fatty acid vesicles^{15,16} and additional justification is provided in the Supplementary Material.
The indirect effect is manifest as a systems property of the model, rather than in any particular function. When a vesicle membrane contains phospholipids (or other surfactant species like hydrophobic peptides), the P_{μ} molecules add a contribution to the surface, increasing the L uptake rate, whereas the L release rate remains unaffected by their presence.
Uptake and release kinetic constants k_{in} and k_{out} are set by two criteria. The first criterion is that pure fatty acid model vesicles (made solely of L), either spherical or deflated, must be in equilibrium when the fatty acid monomer concentration inside and outside the vesicle is the CVC for that amphiphilic compound (e.g. oleic acid). The second criterion is that the model dynamics must reproduce, with lowest RMS error, the experimental time courses reported by Chen et al.^{15} for surface changes in osmotic competition. The second criterion narrows the possible {k_{in}, k_{out}} pairs (see Supplementary Material). For mixed membrane vesicles containing both L and P lipids, we assume that the lipid kinetics equations define what lipid monomer concentration inside and outside the vesicle [L]_{eq} is necessary to keep the mixed membrane vesicle in equilibrium (however, in reality, the CVC of mixed lipid solutions is not a trivial matter^{25}).
For the purpose of lipid competition, has a fixed volume of Ω_{e} litres. Each vesicle has, in principle, a variable internal water volume of litres. This volume value is based on the assumption that water permeates the membrane extremely rapidly and ensures that the interior of each vesicle is isotonic with respect to at all times. However, since in real fatty acid vesicle solutions the concentration of extra compounds like counter ions, pH buffer, etc., is much higher than the concentration of free fatty acid monomers, we can reasonably assume that the aqueous volume of each vesicle is approximately constant at . Thus, vesicle volume is largely determined by the number of buffer molecules a vesicle has trapped inside its internal water pool, with L flux to and from the water pool having marginal osmotic effects.
To summarise, the state of the vesicle system is captured by enumerating the number of lipids in each of the aqueous pools inside the vesicles and in each of the vesicle membranes. The ODE system that describes the time behaviour of the entire vesicle population consists of 2n equations, two for each vesicle:
At the same time, the total number of lipids in the system L_{t} is a conserved quantity set by the initial condition of mixing (see Supplementary Material), always equal to the number of lipid monomers in the environment L_{e}, plus the number of lipids composing the vesicles:
Therefore, L_{e} can be deduced from equation (5) once all L_{c} and L_{μ} have been calculated at time t. Values of model parameters are given in Table 1.
Mean field approximation
In the first instance, before performing any numeric simulations, why should we expect phospholipid fraction and surface growth to be correlated in the vesicle competition model? To answer this question, we can make a mean field approximation. This approach considers a reduced scenario where many details associated to the full model are ignored in order to keep only the logic of the problem (Fig. 3).
The first simplification will be to ignore the internal structure of the vesicles, describing them instead as coarsegrained ‘aggregates’, denoted by pairs , which contain just lipids and phospholipids. This step can be considered justified on the grounds that, at equilibrium, the amount of lipid monomer residing in the vesicle water pools (which typically have tiny volumes, around 1 quintillionth of a litre) is marginal as compared to the lipid composing the vesicle membranes. Since the internal structure or topology of the vesicles is disregarded, it actually amounts to treating them as elongated micelles or flat bilayers.
The second simplification involves reducing the lipid uptake and release equations to their most basic form, independent of membrane tension (u(Φ) = 1) and independent of membrane phospholipid fraction (r(ρ) = 1) respectively. Thus, the ODE system reduces to n simplified equations, where for each aggregate:
Under these conditions, at equilibrium, the molar lipid concentration in the environment [L]_{e} = [L]_{eq} is related to the number of lipids and phospholipids in an aggregate by the following function:
For a fixed number of phospholipids P_{j} > 0, the mapping f:L_{j} → [L]_{eq} can be verified to be onetoone, meaning that each aggregate is in equilibrium at only one specific outside lipid concentration, dependent on the number of lipids L_{j} it contains. Thus, no multiple equilibria of the population are allowed from this type of aggregate dynamics.
Now consider two arbitrarily chosen aggregates i and j in the population of n aggregates, which are competing for lipid. Their ODEs, when written as:
where η = k_{in}/2N_{A}Ω_{e}, are reminiscent of the LotkaVolterra competition equations associated to species sharing and competing for a common set of resources^{26}. If we look for the equilibrium solutions of the previous system, using dL_{i}/dt = dL_{j}/dt = 0, we obtain
which leads to the following proportionality relation at equilibrium:
meaning that the final equilibrium vesicle sizes will be correlated with their respective numbers of phospholipids. Unless P_{i} = P_{j} one of the vesicles will be larger and the second smaller. For each pair (P_{i}, P_{j}) with P_{i} ≠ P_{j} a single solution is found.
When functions u and/or r are not constant, unless they have a trivial form, it is generally not possible to show analytically what shape the correlation between phospholipid fraction and surface growth will take. However, in the Methods section at the end of the paper, we develop a fast numerical way to find the equilibrium configuration of the fullyfledged vesicle population model, with vesicles recovering their internal structure. As compared to numerically integrating the ODE set, the method provides the extra advantages of (i) being faster and thus scaling better for large vesicle populations and (ii) being able to calculate competition ‘tipping points’ (i.e. critical points that mark the transition between growing and shrinking) directly.
In the following Results section, our fast procedure was used to perform accurate vesicle stoichiometry calculations, whilst simulations of the model dynamics were carried out using small populations of vesicles and deterministically integrating the ODE set.
Results
Two competing populations: comparison with experimental results
Figure 4 compares predictions made by our kinetic model against experimentally reported surface growth of vesicles (assessed by a Förster resonance energy transfer assay, FRET) in phospholipiddriven^{17} and osmoticallydriven^{15} competition. Two scenarios are forecast by our model: one where vesicles are spherical when extruded and one where they are deflated by 5% when extruded, as generally observed experimentally^{24,27}. The Supplementary Material details the vesicle population mixing procedure used to initialise our theoretical model, adjusting it to realistic experimental conditions.
Top figures 4a and 4b show the dynamics of surface area change in phospholipiddriven competition. Figure 4a details, in real time, the relative surface area of a tracked (surface area followed by fluorescence probe) population of DOPA:OA (ρ_{0} = 0.1) vesicles, when this population is mixed 1:1 with either pure OA vesicles, similar DOPA:OA (ρ_{0} = 0.1) vesicles or simply buffer. In Fig. 4b, the tracked population is instead pure OA vesicles, which are mixed 1:1 with the same three options outlined above.
Whether starting with initially spherical vesicles, or vesicles deflated by 5%, execution of our lipid kinetics model correctly predicts that when mixed 1:1, DOPA:OA vesicles steal lipid and grow (rising lines, Fig. 4a) at the expense of the pure OA vesicles, which shrink (falling lines, Fig. 4b). In this case, there is also fairly good quantitative agreement with the experimentally observed time courses (RMS error given as Supplementary Table S2). For the other cases, the kinetic model correctly predicts approximately no surface area change (no competition) when similar populations are mixed, or when a population is mixed with buffer.
Middle figures 4c and 4d show phospholipiddriven competition from a different angle: that of vesicle stoichiometry. Stoichiometry explores the final equilibrium size of vesicles in a tracked population, when this population is mixed with a different population containing approximately R times as many vesicles. In this approach, the trend of final equilibrium surface area size versus mixing ratio is explored, rather than the dynamics on the way to equilibrium. Figure 4c details final surface area of a tracked population of DOPA:OA (ρ_{0} = 0.1) vesicles, when this population is mixed 1:R with a population of pure OA vesicles. Figure 4d details the opposite scenario, whereby the tracked population is OA vesicles, mixed 1:R with DOPA:OA vesicles. The R = 1 cases in Figs 4c and 4d correspond to the surface sizes reached in the limit of time in Figs 4a and 4b, respectively.
Calculating competition equilibrium by means of the fast computation approach outlined in the Methods section, we were able to verify that our model exhibits continual growth of DOPA:OA (ρ_{0} = 0.1) vesicles as more OA vesicles are added (Fig. 4c). If vesicles started at 5% deflation, the model matched the experimental data points even more closely. In the opposite scenario, we also verified that the model shows the same distinctive plateau in the shrinkage of pure OA vesicles as more DOPA:OA (ρ_{0} = 0.1) vesicles are added (Fig. 4d). In the case of the latter figure, notably the indirect effect alone is sufficient to reproduce experimental results.
Importantly, the general outcome of phospholipiddriven competition in our model is for DOPA:OA mixed vesicles to uptake lipid, grow in surface and to finish at high Φ > 1 values (excess surface, flaccid), whereas pure OA vesicles lose lipid, suffer reduced surface and finish at Φ < 1 values (osmotically tense, spherical). This is observed experimentally and indeed provides the basis for the conjecture that phospholipidladen vesicles are more likely to divide spontaneously when gentle external shearing forces are applied^{17}.
Moving to osmoticallydriven competition, Fig. 4e shows simulation of a swelled population of vesicles competing with an initially isotonic (nonswelled) population. Simulation outcomes match quite well the experimental bestfit time courses, in particular for the growth of the swelled vesicles (less accurately for the shrinkage of the nonswelled vesicles). In any case, it must be noted that the original experimental data (yellow data points) has considerable variance. Then, Fig. 4f shows that the kinetics model qualitatively reproduces the stoichiometric observation whereby adding more swelled vesicles to a population of initially nonswelled vesicles will cause the shrinkage of the nonswelled vesicles to plateau, rather than to continue (note the logarithmic scale of Fig. 4f). Again, model outcomes are improved if vesicles start at 5% deflation.
The general outcome of osmoticallydriven competition in our model is for all vesicles to finish with different surface sizes (as for phospholipiddriven competition), but now, all vesicles also share the same Φ < 1 value, indicating equal osmotic stress. This residual osmotic stress is also observed experimentally and stands as the main criticism of the osmoticallydriven competition scenario. In order to divide, swelled vesicles would have to overcome a stronger energetic barrier, changing their stressed membrane state into one ready for fission, making this an improbable route to spontaneous vesicle reproduction^{18}.
Our kinetic model can also be used to make predictions or to find competition ‘tipping points’ in the more general scenario where completely heterogeneous populations of phospholipidladen and/or osmotically swollen vesicles compete for lipid (Figs. 5 and 6), even if some of these experiments have not been realised in the lab yet.
Competition tipping points in diverse populations
Figure 5a shows that within a hypothetical population of model phospholipidladen vesicles, where each vesicle has a randomly assigned phospholipid fraction in the membrane between 0 and 100%, the critical DOPA fraction needed for growth (tipping point), in this case, is just over 58%.
Figure 5b compares different heterogeneous populations competing for phospholipid and reveals an important observation: competition is always context dependent. That is to say, a certain amount of membrane phospholipid does not guarantee a certain final surface area. Rather, final surface depends on the boundary conditions of the competition event (that is, the parameters influencing the solution of equation (15) in the later Methods section), which includes the number and composition of competitor vesicles present. For example, population (i) in Fig. 5b has vesicles with low DOPA fraction as compared to vesicles in population (iv), yet in some cases the vesicles in the former population have larger final surface growth than vesicles in the latter. This concurs with the experimental observation that even small differences in phospholipid content should drive growth^{17}.
The dotted black lines in Figs 5a and 5b are the same competition events run when the direct effect is present and maximally enabled (d = 1). The extent to which the direct effect affects vesicle growth must be made on a case by case basis, as it depends on the specifics of the competition event. For example, the direct effect has marginal influence on vesicle growth trends in the population shown in Fig. 5b (iii), but is more relevant in population (ii). The Supplementary material contains a recalculation of both Fig. 5a and 5b, if we further take into account the realistic constraint that vesicles will burst when osmotic tension exceeds a critical limit (Φ < 0.7).
Figure 5c shows that in a heterogeneous population where pure OA model vesicles are swelled to differing extents, vesicles with low initial Φ values take lipid from those with higher (less swelled) Φ values, with the tipping point between growing and shrinking at . As a last remark, orange crosses marked on Figs 5a and 5c show that full deterministic simulations of the model (run all the way to equilibrium) agree with and thus validate the ‘Fast Computation of Competition Equilibrium’ procedure outlined in the Methods section.
Theoretical predictions beyond current experimental results
Finally, we were able to explore more widely some of the parameter space for phospholipiddriven and osmoticallydriven competition, using our model to make some predictions. Figure 6a shows the stoichiometry results of phospholipiddriven competition in this wider context. A population of DOPA:OA (ρ_{0} = 0.1) vesicles is mixed with a second population, but the phospholipid content of the second population, as well as the mixing ratio R, are varied. Taking a slice through the surface labelled ‘pop1’ when shows the result reported as the solid red line in Fig. 4c. Figure 6b explores the stoichiometry of osmoticallydriven competition in a similar way to phospholipiddriven competition. A fixed population of swelled vesicles is mixed with a second population, where the degree of swelling in the second population, as well as the mixing ratio R are varied. To conclude these predictions, Fig. 6c shows the effects of osmoticallydriven versus phospholipiddriven competition, still a completely unreported scenario in the experimental literature, whereby a population of swelled pure oleate vesicles competes for lipid with a population of DOPA:OA vesicles. The swelled oleate vesicles are able to steal lipid from the DOPA:OA vesicles, when the former have a high degree of swelling and the latter have a low DOPA fraction; otherwise, the DOPA:OA vesicles prosper in the competition.
Vesicle bursting is an important consideration in Fig. 6. Competition predictions in Fig. 6a are only strictly valid when the population 1 surface is above the red box lines. Below these lines, vesicles in population 1 have excessive osmotic pressure (Φ < 0.7) and would likely burst, altering competition outcomes for population 2. Likewise, the population 2 surface in Fig. 6c is only drawn for values where oleate vesicles in that population have final Φ > 0.7. Outside the extent of the population 2 surface, competition outcomes for population 1 should be treated with caution, as not all population 2 vesicles will be intact.
Discussion
In this work we have presented a theoretical model of the transfer kinetics of single chain fatty acids between competing vesicles. We have shown that data coming from controlled laboratory experiments on phospholipiddriven competition and osmoticallydriven competition can be reproduced fairly well by a set of physicallybased rate equations describing the uptake and release of fatty acids for each vesicle. Furthermore, we have been able to predict the outcome of several yettobeperformed experiments. Thus, it is time to recapitulate, considering possible limitations of our approach, clarifying several points that remain open and giving a more general perspective on the problem addressed.
The main assumption we made when modelling phospholipiddriven competition is that the phospholipids are not released by vesicle membranes at the timescale of fatty acid transfer between the supramolecular structure (i.e., the closed membrane bilayer) and the aqueous solution (both inwards and outwards). This assumption is well founded on experimental evidence^{14,15} and a previous theoretical model^{16}. Likewise, the assumption we make with osmoticallydriven competition is that the extra buffer present inside the vesicles, swelling them, permeates very slowly through the bilayer membrane. In reality, the offrate of a lipid molecule from a bilayer is inversely correlated with the number of carbon atoms in the acyl chain of the lipid concerned and phospholipids do have a small nonzero transfer rate (with a half time from hours to days^{16,17}). If the much slower phospholipid transfer was included in our model, the equilibrium reached in the limit of time would always be that of a completely homogeneous population. This is because the P phospholipid would redistribute amongst the vesicles until all were equilibrated with the same phospholipid monomer concentration in solution [P]_{eq}, which is trivially when all vesicles have the same fraction of membrane phospholipids P_{μ}. With no remaining asymmetries in P_{μ} fraction to drive competition, all vesicles would finish with the same lipid composition and same surface size. The initial appearance and eventual disappearance of competition would thus follow the same type of dynamics as those experimentally reported by Budin & Szostak^{17} (Fig. 3D therein) for nervonic acid, which redistributes between vesicles. If vesicles contained a metabolism which synthesised phospholipid, then lasting P_{μ} asymmetries between vesicles could conceivably be maintained as steady states, despite a continuous process of P exchange. However, in this work we took the route of not explicitly modelling phospholipid synthesis, to reduce the competition scenario to a materiallyclosed system which subsequently settles to equilibrium. Under this condition, analysis is easier to perform. In summary, the results of this study can be interpreted as reflecting the competition advantage bestowed upon a vesicle by a membrane phospholipid fraction given that this fraction is somehow maintained as constant.
The next point that deserves discussion is the the causative role of the direct effect in the phospholipiddriven competition simulations performed with our model. Our choice for function r means that a DOPA:OA vesicle with 10 mol% DOPA fraction will have fatty acids leaving the bilayer at reduced rate r(0.1) × k_{out} = 0.9k_{out} when the direct effect is maximal (d = 1). It could be argued that other function choices for r could reduce fatty acid offrate even further for the same DOPA fraction. However, a large direct effect is not needed to best fit model outcomes with experimental outcomes and would actually make the fit worse. Examination of Figs 4b and 4d in fact shows that having only the indirect effect provides the best fit to experimental outcomes (quantified in Supplementary Table S2). On the other hand, the dynamics and stoichiometry outcomes of Figs 4a and 4c respectively are only improved when there is a small reduction of k_{out}: a small direct effect of around 0 < d < 1 for Fig. 4a and with d a little larger than 1 for Fig. 4c. With the maximal level of direct effect provided by our function r, Supplementary Table S3 shows that the direct effect only accounts for around 20% of the total vesicle surface growth. Therefore, we should conclude that in our kinetic treatment of vesicle competition, the indirect effect is the main mechanism driving vesicle growth dynamics.
One curiosity in the results (both in vitro and in silico) is how DOPA:OA (ρ_{0} = 0.1) vesicles grow continually as more OA vesicles are added (Fig. 4c). This is unintuitive, since the growth of the DOPA:OA vesicles should imply a dilution of their phospholipid content, which would seemingly reduce the indirect and direct effects, thus giving a negative feedback to eventually curb the DOPA:OA growth profile. The reason why our model reproduces this continuous growth result has to do with the mathematics underlying the kinetic modelling. In the limit of infinite L_{μ} lipids in the membranes of our model DOPA:OA vesicles, the inside/outside lipid concentration required to sustain them at equilibrium (given by function f, as defined by equation (13) later in the Methods section) tends to, but crucially never actually reaches, the CVC of pure oleic acid:
This is true, even if a model DOPA:OA vesicle contains just one single phospholipid in the membrane. Now, as more OA vesicles are mixed with the DOPA:OA vesicles, the population becomes increasingly dominated by OA vesicles and the lipid monomer concentration in the environment subsequently rises toward . As this happens, equation (10) implies that the DOPA:OA vesicles will be absorbing more and more L lipids, in order to grow to a size that is in equilibrium with the external lipid monomer concentration. The growth of our model DOPA:OA vesicles is thus halted only by the number of lipids in the system being limited to L_{t} and not by dilution of the membrane phospholipid fraction. In our kinetics model, continuous growth happens with or without the direct effect present.
A final point worth highlighting is that when the lipid uptake function u given in equation (2) is not conditional, as we assumed, but simply
for all membrane states (which denotes that even flaccid vesicles have differential rates of lipid uptake), then, quite interestingly, the continuous DOPA:OA growth effect cannot be reproduced. If we use definition (11) for function u in function f (13) and call the new function f_{nc}, it can be shown that
meaning that the DOPA:OA vesicles do not show the same continued growth as the lipid monomer concentration in the environment rises toward . Rather, the DOPA:OA have much slower growth and they even have a finite stable size when the outside lipid monomer concentration is exactly . Thus, to best reproduce experimental outcomes, a crucial part of our lipid uptake kinetics was to accelerate lipid uptake only in osmotically stressed vesicle states, not in flaccid ones. This is a new addition to our general kinetic model, introduced in this paper.
This work is a step forward in the development of semirealistic, coarsegrained descriptions of phenomena that, in reality, are extremely complex. Selfassembly processes involving heterogeneous component mixtures and the formation of dynamic supramolecular structures that could hypothetically lead to biologically relevant forms of material organization, like protocells^{28,29}, constitute a tremendous challenge, indeed, both for experimental and theoretical ‘systems chemistry’ research^{30} and for synthetic biology^{7,31}. In particular, the connection between basic metabolic reaction networks and membrane dynamics (including stationary growth and division cycles^{32}) needs to be explored much more extensively, since it is one of the key aspects to establish a plausible route from physics and chemistry towards biological phenomenology.
Methods
Fast computation of competition equilibrium
Here we provide a general numerical approach to solving the equilibrium configuration of a possibly heterogeneous population of vesicles competing for a limited supply of lipid. These vesicles may be osmotically swelled, laden with phospholipid, or a mixture of both and can be arbitrary in number. The method allows the lipid uptake and release functions u and r to take arbitrary forms, subject to some requirements detailed below.
We start by defining a function f:L_{μ} → [L]_{eq}, like equation (7), which gives the inside/outside lipid monomer concentration [L]_{eq} necessary to maintain a particular vesicle at equilibrium, given that this vesicle has a specific number of lipids/phospholipids in the membrane and a specific volume:
The inverse of this function yields useful information: it is the mapping of [L]_{eq} to the number of lipids which must exist in the membrane of a particular vesicle, in order for that vesicle to be at equilibrium.
However, due to the difficulty in isolating L_{μ} from the potentially nonlinear functions u and r, in most cases the inverse mapping is not possible to write in closed form. Nevertheless, if uptake and release functions u and r make function f (i) onetoone, (ii) onto and (iii) continuous, then it follows that the inverse mapping is a function f^{−1}, which can be numerically calculated for vesicle by using f and binary searching for an L_{μ} which satisfies:
using appropriate search bounds (normally: , ).
Crucially, having a means to calculate f^{−1} gives a way of determining the total number of lipids existing in all equilibrated vesicle membranes, given that the inside/outside lipid monomer concentration in the heterogeneous vesicle mixture is [L]_{eq}. For each [L]_{eq}, we know that each vesicle has a unique number of membrane lipids L_{μ}, because f^{−1} is itself onetoone. This means that a certain [L]_{eq} can only admit one single equilibrium configuration of vesicles, not multiple equilibrium configurations and this lack of ambiguity is a desirable property for the method.
The lipid monomer concentration [L]* inside/outside all vesicles in this single equilibrium configuration can be found by making use of the lipid conservation principle in equation (5):
That is, at [L]*, the lipid making up the membranes of all equilibrated vesicles, plus the lipid monomer inside and outside the vesicles is equal to the total lipid in the system L_{t} set by the initial condition. Expression (15) can also be solved by binary search of [L]_{eq} between appropriate bounds, normally over all vesicles, [L]_{eq}^{max} = min(f(L_{μ} = L_{t})) over all vesicles. Finally, knowing [L]* allows to fully reconstruct the final sizes of all vesicles at equilibrium by substituting [L]_{eq} = [L]* into equation (14) for each vesicle .
In the equilibrated population, some vesicles will have grown larger in surface area at the expense of others which will have shrunk. When a population of vesicles has competed for lipid via phospholipiddriven competition, the ‘tipping point’ is the critical number of membrane phospholipids separating those vesicles which have lost lipid from those which have gained lipid and is found by:
where , again solvable by binary searching, this time in the range , (from a pure lipid membrane to a pure phospholipid membrane). Expression (16) amounts to asking how many phospholipids a hypothetical vesicle would require in order not to grow in surface area when the lipid monomer concentration has stabilised at [L]*. Likewise, the number of phospholipids required to achieve any arbitrary surface area growth can be found by setting S_{μ} to the value desired.
The critical phospholipid number can be stated more usefully as the critical phospholipid molecular fraction
a vesicle has in the initial condition, a time when all vesicles have a surface of . For osmoticallydriven competition, the critical volume separating shrinking vesicles from growing vesicles is found by searching expression (16) for vesicle volume instead:
where . This may be alternatively stated as the critical Φ in the initial condition:
If no sign change results when evaluating the functions given in equations (14–18) at the upper and lower search bounds, then the respective equation cannot be solved by this numerical bisection approach. Otherwise typically 30 iterations of binary search were used to converge to an accurate answer.
References
Smith, J. M. & Szathmáry, E. The Major Transitions in Evolution (Oxford University Press, 1995).
Luisi, P. L. The Emergence of Life: From Chemical Origins to Synthetic Biology (Cambridge University Press, 2006).
Varela, F. G., Maturana, H. R. & Uribe, R. Autopoiesis: the organization of living systems, its characterization and a model. Biosystems 5, 187–196 (1974).
Gánti, T. Organization of chemical reactions into dividing and metabolising units: the chemotons. Biosystems 7, 15–21 (1975).
Dyson, F. Origins of Life (Cambridge University Press, 1985).
Segré, D., BenEli, D., Deamer, D. W. & Lancet, D. The lipid world. Orig Life Evol Biosph 31, 119–145 (2001).
Solé, R. V., Munteanu, A., RodriguezCaso, C. & Macía, J. Synthetic protocell biology: from reproduction to computation. Philos Trans R Soc Lond B Biol Sci 362, 1727–1739 (2007).
Macía, J. & Solé, R. V. Synthetic turing protocells: vesicle selfreproduction through symmetrybreaking instabilities. Philos Trans R Soc Lond B Biol Sci 362, 1821–1829 (2007).
Mavelli, F. & RuizMirazo, K. Stochastic simulations of minimal selfreproducing cellular systems. Philos Trans R Soc Lond B Biol Sci 362, 1789–1802 (2007).
RuizMirazo, K., Piedrafita, G., Ciriaco, F. & Mavelli, F. Stochastic Simulations of MixedLipid Compartments: From SelfAssembling Vesicles to SelfProducing Protocells. vol. 696 of Adv Exp Med Biol, 689–696 (Springer, 2011).
Mavelli, F. Stochastic simulations of minimal cells: the ribocell model. BMC Bioinformatics 13, S10 (2012).
Munteanu, A., Attolini, C. S.O., Rasmussen, S., Ziock, H. & Solé, R. V. Generic darwinian selection in catalytic protocell assemblies. Philos Trans R Soc Lond B Biol Sci 362, 1847–1855 (2007).
Mavelli, F. & Stano, P. Kinetic models for autopoietic chemical systems: the role of fluctuations in a homeostatic regime. Phys Biol 7, 16010 (2010).
Cheng, Z. & Luisi, P. L. Coexistence and mutual competition of vesicles with different size distributions. J Phys Chem B 107, 10940–10945 (2003).
Chen, I. A., Roberts, R. & Szostak, J. W. The emergence of competition between model protocells. Science 305, 1474–1476 (2004).
Mavelli, F. & RuizMirazo, K. Environment: a computational platform to stochastically simulate reacting and selfreproducing lipid compartments. Phys Biol 7, 036002 (2010).
Budin, I. & Szostak, J. Physical effects underlying the transition from primitive to modern cell membranes. Proc Natl Acad Sci U S A 108, 5249–5254 (2011).
Adamala, K. & Szostak, J. W. Competition between model protocells driven by an encapsulated catalyst. Nat Chem 5, 495–501 (2013).
Erdi, P. & Tóth, J. Mathematical models of chemical reactions: Theory and applications of deterministic and stochastic models (Manchester University Press, 1989).
Ullah, M. & Wolkenhauer, O. Stochastic Approaches for Systems Biology (Springer, 2011).
Simard, J. R., Pillai, B. K. & Hamilton, J. A. Fatty acid flipflop in a model membrane is faster than desorption into the aqueous phase. Biochemistry 47, 9081–9089 (2008).
Kamp, F. & Hamilton, J. A. ph gradients across phospholipid membranes caused by fast flipflop of unionized fatty acids. Proc Natl Acad Sci U S A 89, 11367–11370 (1992).
Chen, I. A. & Szostak, J. W. Membrane growth can generate a transmembrane ph gradient in fatty acid vesicles. Proc Natl Acad Sci U S A 101, 7965–7970 (2004).
Mally, M., Peterlin, P. & Svetina, S. Partitioning of oleic acid into phosphatidylcholine membranes is amplified by strain. J Phys Chem B 117, 12086–12094 (2013).
Cape, J. L., Monnard, P.A. & Boncella, J. M. Prebiotically relevant mixed fatty acid vesicles support anionic solute encapsulation and photochemically catalyzed transmembrane charge transport. Chem Sci 2, 661–671 (2011).
Lotka, A. J. Elements of Physical Biology (Williams and Wilkins, Baltimore, 1925).
Mui, B. L.S., Cullis, P. R., Evans, E. A. & Madden, T. D. Osmotic properties of large unilamellar vesicles prepared by extrusion. Biophys J 64, 443–453 (1993).
Mouritsen, O. G. Life  As a Matter of Fat: The Emerging Science of Lipidomics (SpringerVerlag, 2005).
Rasmussen, S. et al. (eds.) Protocells: Bridging NonLiving and Living Matter (MIT Press, 2009).
RuizMirazo, K., Briones, C. & de la Escosura, A. Prebiotic systems chemistry: New perspectives for the origins of life. Chem Rev 114, 285–366 (2014).
Solé, R. V. Evolution and selfassembly of protocells. Int J Biochem Cell Biol 41, 274–284 (2009).
Mavelli, F. & RuizMirazo, K. Theoretical conditions for the stationary reproduction of model protocells. Integr Biol 5, 324–341 (2013).
Lonchin, S., Luisi, P. L., Walde, P. & Robinson, B. H. A matrix effect in mixed phospholipid/fatty acid vesicle formation. J Phys Chem B 103, 10910–10916 (1999).
Acknowledgements
R.S. and B.S.E. acknowledge support from the Botin Foundation and by the Santa Fe Institute. K.R.M. acknowledges support from the Basque Government (Grant IT 59013), Spanish Ministry of Science (MINECO Grant FFI201125665), COST Action CM 1304 (Emergence and Evolution of Complex Chemical Systems). F.M. acknowledges support from MIUR (PRIN 2010/11 2010BJ23MN_003). We thank Itay Budin and Irene Chen for kindly providing original experimental data for Figure 4 and the Group of Dynamical Systems (Department of Applied Mathematics and Analysis) from Universitat de Barcelona for providing us with the RungeKuttaFehlberg algorithm used for numerical integration of the model.
Author information
Affiliations
Contributions
All authors helped conceive the model. B.S.E. and R.S. analysed the model and B.S.E. conducted the numeric simulations. B.S.E., K.R.M. and R.S. wrote the paper. F.M. provided crucial feedback on several aspects.
Ethics declarations
Competing interests
The authors declare no competing financial interests.
Electronic supplementary material
Supplementary Information
Supplementary Material for: Modelling Lipid Competition Dynamics in Heterogeneous Protocell Populations
Rights and permissions
This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder in order to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/
About this article
Cite this article
ShirtEdiss, B., RuizMirazo, K., Mavelli, F. et al. Modelling Lipid Competition Dynamics in Heterogeneous Protocell Populations. Sci Rep 4, 5675 (2014). https://doi.org/10.1038/srep05675
Received:
Accepted:
Published:
Further reading

Prebiological Membranes and Their Role in the Emergence of Early Cellular Life
The Journal of Membrane Biology (2020)

Is defining life pointless? Operational definitions at the frontiers of biology
Synthese (2018)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.