Modelling Lipid Competition Dynamics in Heterogeneous Protocell Populations

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 cell-level 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 coarse-grain 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.


I. INTRODUCTION
A fundamental problem in biology involves the origins of an innovation that allowed the development of organisms in our biosphere, beyond complex chemical reaction networks: the emergence of cells (Luisi, 2006;Smith and Szathmáry, 1995).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.However, although it is often assumed that there was a transition from some type of 'lessorganised' prebiotic chemistry (surely including catalytic cycles) to a cell-based living chemistry, little is yet known concerning 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 (Dyson, 1985;Gánti, 1975;Macía and Solé, 2007;Mavelli and Ruiz-Mirazo, 2007;Ruiz-Mirazo et al., 2011;Segré et al., 2001;Solé et al., 2007;Varela et al., 1974).In particular, by considering the coupling of a tem-plate carrying information with vesicle replication and metabolism, it has been shown that Darwinian selection is the expected outcome of competition in a protocellular world (Munteanu et al., 2007).However, an early pre-Darwinian stage in the development of biological organisms was likely to be dominated by supramolecular systems disconnected from information, closer to elementary forms of metabolism and strongly constrained by the molecular diversity of the available chemical repertoire.What type of competition and cooperation processes were at work in the chemical world leading to the emergence of early protocells?Here, in that context, processes able to favour asymmetries in the chemical composition of vesicles should be expected to play a relevant role in defining the conditions under which protocellular assemblies could thrive.
Recent laboratory experiments have actually shown how simple physical changes made to the lipid membrane of vesicles can drive competition between those vesicles when the supply of lipid is limited.First, (Chen et al., 2004) reported competitive dynamics in a population of 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 lipids more quickly.More recent experimental work has turned attention to other possible selective advantages of protocells, such as phospholipid- (Budin and Szostak, 2011) and peptide- (Adamala and Szostak, 2013) driven competition amongst vesicles.Instead of membrane tension, the main factor for competition here is a different type of molecule inserted in the lipid bilayer, which changes its physical properties.In Two mechanisms of phospholipid-driven growth.A Indirect effect, whereby the presence of phospholipid in a vesicle membrane drives growth simply through a geometric asymmetry: only the lipid section of the bilayer (grey) is able to release lipid (orange arrows) whereas the whole of the bilayer surface (made of lipids and phospholipids) is able to absorb lipid monomer (green arrows).Phospholipid fraction is pictured as one continuous block to highlight the principle only.The indirect effect can be created also by nonlipid surfactant molecules (e.g.peptides) residing long enough in the membrane to increase surface absorption area.B Direct effect, whereby the acyl tails of the phospholipids have high affinity for packing closer to each other and increasing bilayer order, thus making the exit of the simple lipids more difficult.The direct effect is specific to the molecular structure of phospholipids.In both cases, growth eventually stops when the phospholipid fraction in the membrane becomes diluted.
the first case, which will be the main focus of this work, fatty acid vesicles endowed with a membrane fraction of phospholipid are observed to steal lipid molecules from phospholipid-deficient neighbours, who shrink, whilst the former grow and keep their potential for division.Two basic physical mechanisms have been postulated to underlie phospholipid-driven growth, as explained in Fig. 1 under the terms indirect and direct effects.
With the aim to complement experimental results, and in an attempt to better formalise and investigate competition processes at play, in this paper we develop a mathe-matical model of a competing population of vesicles.The model is based at the level of lipid kinetics, following the approach of (Mavelli and Ruiz-Mirazo, 2010).A vesicle in the population absorbs and releases lipid 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 swelling).Using physically realistic parameters (i.e.lipid molecule sizes, vesicle aggregation numbers and CVC concentrations -see Table I) we are able to qualitatively and often quantitatively reproduce experimental data for phospholipid-driven and osmotically driven competition.
The paper is organised as follows.The Methods section introduces the kinetic model.A mean-field analysis is performed to give insight into why we should expect phospholipid-driven competition to result from a basic version of the model kinetics, followed by the description of a fast numerical method for solving the final equilibrium state of the full model.Then, the vesicle mixing procedure is specified, in order to be able to interface the model with experimental observations.The Results section summarises how well the kinetic model is able to reproduce experimental results and observations, including also some predictions for still untested protocell competition scenarios.Finally, in the Discussion section, we consider several possible limitations of our approach and conclude the study.

Author Summary
Synthetic protocell biology is bringing forward exciting experimental results that allow us to conceive in more realistic terms how the first living organisms could have emerged and started a process of Darwinian evolution.A remarkable finding has been the capacity of lipid vesicle populations to undergo competition and selection processes without the need of nucleotide replication mechanisms.This opens a completely new research avenue to explore and characterize pre-Darwinian modes of evolution leading to the first protocellular systems.In this work, we develop a mathematical model of vesicle competition to complement ongoing experimental efforts and to also provide a reliable way to investigate scenarios or conditions that are difficult to survey in the wet lab.Our model, which is based at the level of lipid kinetics, is demonstrated to reproduce diverse reported results and is helpful in providing new insights about the molecular mechanisms underlying protocell growth and competition dynamics.

A. Theoretical Model of Vesicle Competition
The competition model involves a set of n vesicles embedded in a finite volume environment E defined by a triple (Ω e , L e , B e ).Each competing vesicle consists of a unilamellar membrane of two different lipid types: a fixed number of phospholipids P µ (e.g.di-oleoylphosphatidic acid, DOPA) and a variable number of simple fatty acid lipids L µ (e.g.oleic acid, OA).The L lipids in the bilayer continuously exchange with the vesicle internal water pool (also considered a well-mixed chemical domain), and E, whereas the P phospholipids are considered approximately stationary due to their comparatively slow exchange rate (a reasonable assumption, given the very low CVC values of standard phospholipids compared to other species in water).The internal water pool of each vesicle hosts L c lipid monomers and also B c buffer species, which cannot permeate the bilayer but provide osmotic stability.These buffer species are also present in E with constant number B e .
Vesicles compete with each other by consequence of uptaking/releasing simple lipids L from/to E, 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 non-equilibrium 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 vesicle-vesicle interactions, and assume a well-mixed set of vesicles (Fig. 2).
More precisely, each vesicle V j is considered to release lipids to both aqueous phases (at each side of the bilayer) at the equal rate of λ out j = k out L j µ r(ρ j ), and absorb lipids from each phase at rate λ in j = k in S j µ [L]u(Φ j ), where [L] is the molar concentration of lipid monomer in the respective phase.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 [L] j c = [L] e = [L] * at equilibrium.Flip-flop of the simple lipid L between membrane leaflets is considered very fast with respect to its uptake and release rates, and thus the bilayer is modelled as a single oily phase.
The total number of lipids in the system L t is a conserved quantity set by the initial condition of mixing, always equal to the number of lipid monomers in the environment L e , plus the number of lipids composing the vesicles.Therefore, at all times: The state of the system is captured by enumerating the number of lipids in each of the aqueous pools inside the vesicles, and each of the vesicle membranes.The ODE system consists of 2n equations, two for each vesicle: and L e can be deduced from constraint (1), once all L c and L µ have been calculated at time t.
Explaining the choice of lipid L release kinetics, each lipid in a pure L membrane is considered to have a uniform probability per unit time k out of disassociating from the membrane, and function r modifies this probability based on the current molecular fraction of phospholipid in the membrane ρ = Pµ Pµ+Lµ .In order to account for the direct effect, we define function 0 ≤ r(ρ) ≤ 1 to be monotonically decreasing with increasing ρ, meaning that increasing phospholipid fraction generally decreases bilayer fluidity, slowing down the rate of L release from the membrane (Budin and Szostak, 2011).In a first approximation, r is 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 geometric factor Φ = S µ / 3 √ 36πΩ 2 (where Ω is vesicle aqueous volume, in the same units as S µ ).We define a conditional function to denote that lipid uptake is only increased when the the bilayer is stressed (Φ < 1).Flaccid vesicles do not have extra enhancement of lipid uptake rate.Rationale for this function stems from (Mavelli and Ruiz-Mirazo, 2010), to account for osmotically-driven competition between vesicles (Chen et al., 2004).
To clarify some final assumptions, the vesicle surface area, referred to as S µ = 1 2 (L µ α L + P µ α P ), is the waterexposed area of the inner bilayer leaflet, or alternatively the water-exposed area of the outer bilayer leaflet of the 2 Kinetic model of vesicle competition.A Our model approach considers as a starting point a population of vesicles (of generally heterogeneous sizes and membrane compositions) in a well-mixed environment.B Each vesicle has a membrane composed of simple single chain lipids L, e.g.oleic acid (OA), and C sometimes more complex double chain phospholipids P , e.g.di-oleoyl-phosphatidic acid (DOPA).D outlines the kinetic interactions between vesicles.Here two vesicles are displayed.
Vesicle 1, on the left hand side, has a mixed membrane with approximately 10 mol % phospholipids P (black) and the remainder single chain lipids L (grey).Vesicle 2 consists purely of simple lipids L. In the ensuing competition, phospholipid-laden vesicle 1 will grow at the expense of vesicle 2, which will shrink.
vesicle.Membrane thickness is therefore considered negligible.Uptake and release kinetic constants k in and k out are set taking into account that spherical vesicles made purely of L should be in equilibrium when the lipid monomer concentration inside and outside the vesicle is the experimentally observed CVC concentration for that amphiphilic compound (e.g.oleic acid), and also considering that L uptake is orders of magnitude faster than L release.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 (Cape et al., 2011)).
For the purpose of lipid competition, E has a fixed volume of Ω e litres, and each vesicle V j has, in principle, a variable volume internal water pool of Ω j = Ω e (L j c + B j c )/(L e + B e ) litres.This condition ensures that, at all times, the interior of each vesicle is isotonic with respect to E. However, we make the simplifying assumption in this work that vesicles exist in a solution with a comparatively high buffer concentration.Thus, each vesicle has an approximately constant aqueous volume Ω j ≈ Ω e (B j c /B e ) largely determined by the number of buffer molecules it has trapped inside the internal water pool, with L flux to and from the water pool having marginal osmotic effects.Model parameters are given in Table I.

B. Mean Field Approximation
With the goal to gain intuition about why one should expect phospholipid fraction and surface growth to be correlated in the vesicle population model described, we can make a mean field approximation.This approach considers a reduced scenario where many details associ-... ...
3 Meanfield model of vesicle population dynamics.Considering vesicles as simplified aggregates permits some analytical treatment.
ated 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 coarse-grained 'aggregates', denoted by pairs V j = (L j , P j ), 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 one-to-one, 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 Lotka-Volterra competition equations associated to species sharing and competing for a common set of resources (Lotka, 1925).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: This result immediately tells us that, for a given fraction P i /P j the relative sizes of the two chosen vesicles are correlated.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 next section we develop a fast numerical way to find the equilibrium configuration of the fully-fledged 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.

C. Fast Computation of Competition Equilibrium
In this section 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 (7), which gives the inside/outside lipid monomer concentration [L] eq necessary to maintain a particular vesicle V j 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 non-linear 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 both (i) one-to-one and onto and (ii) continuous, then it follows that the inverse mapping is a function f −1 , which can be numerically calculated for vesicle V j 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 one-to-one.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 (1): 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 (12) can also be solved by binary search of [L] eq between appropriate bounds, normally [L] (min) eq = max(f (L µ = 0)) 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 (11) for each vesicle V j .
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 phospholipid-driven competition, the 'tipping point' is the critical number of membrane phospholipids P crit µ separating those vesicles which have lost lipid from those which have gained lipid, and is found by: where S µ = S 0 µ , again solvable by binary searching, this time in the range 0 ≤ P µ ≤ 2S 0 µ /α P , (from a pure lipid membrane to a pure phospholipid membrane).Expression (13) 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 S 0 µ .For osmotically-driven competition, the critical volume separating shrinking vesicles from growing vesicles is found by searching (13) for vesicle volume instead: where S µ = S 0 µ .This may be alternatively stated as the critical Φ in the initial condition: If no sign change results when evaluating the functions (11 -15) 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 1 .
1 Model vesicles described in this work either contain just fatty acid membranes, or fatty acid mixed with one other type of phospholipid.Another avenue (not explored here) would be to vary the type of phospholipid from vesicle to vesicle.In this case, competition equilibrium can still be calculated by adding two

D. Modelling Vesicle Mixing
In order to interface our theoretical model with experimentally reported results, it is necessary to define an acceptable procedure for mixing two competing vesicle populations.The basic mixing procedure outlined in this section establishes the boundary conditions for competition, which are (i) the number of vesicles present, (ii) their respective compositions, (iii) the environment volume Ω e and (iv) the total amount of lipid L in the system (L t ).Below, turning more specific to match experimental scenarios, simple lipids L are considered to be OA, and phospholipids P are considered to be DOPA.

Competition Volume
The equilibrium finding method outlined in the 'Fast Computation of Competition Equilibrium' section above requires summing over a finite number of vesicles.Likewise, dynamic simulations of the model require a finite ODE set.However, vesicle populations in a real laboratory experiment will typically have millions of vesicles competing for lipid.In our modelling approach, it is therefore necessary to consider a small volume 'patch' of each of the solutions being mixed.Each patch volume is large enough to contain enough vesicles so as to be representative of the vesicle density in the solution it pertains to, but no so many vesicles that numerical solution becomes infeasibly slow.
A patch volume Ω p = Ω stoi litres (Table I) was utilised for stoichiometric calculations using the equilibrium finding method, which translates into around 2000 vesicles being involved in 1:1 mixing.Full dynamic simulation of the model with the Gillespie Direct SSA algorithm forced a yet smaller patch volume Ω p = Ω dyn litres to be used (to give τ jumps which were not too small), translating into around 40 vesicles being involved in 1:1 mixing.

Mixing for Phospholipid-Driven Competition
In order to mix a fixed population of DOPA:OA vesicles which have molecular fraction ρ of DOPA in their membranes, in ratio R with a variable population of pure OA vesicles, we assume the following basic steps.
Firstly, a suspension of OA lipid monomers in concentration RC 0 molar (assuming RC 0 CVC for oleic acid) is extruded (possibly multiple times) through 100nm diameter pores.This leads to a more homogeneous population of 100nm diameter pure OA unilamellar vesicles.more arguments to function f : the first would detail the head area for the phospholipid in vesicle V j ; the second would describe how this phospholipid changes bilayer fluidity and alters the simple lipid L release rate (e.g.parameter d to function r could be supplied).
Each vesicle is assumed spherical (Φ = 1) with aqueous volume Ω 0 .The molar concentration of OA vesicles in the extruded suspension is approximately where N OA is called the 'aggregation number', equal to the total number of lipids forming a vesicle bilayer (in this case, just OA lipids).The lipid monomer concentration in the aqueous solution inside/outside the vesicles is [L] OA eq , the CVC value, maintaining them at equilibrium.Secondly, a mixed suspension containing both OA lipid monomers (in molar concentration C 0 ) and DOPA phospholipids (in molar concentration gC 0 , where g = ρ 1−ρ ) is extruded through 100nm diameter pores.This, similarly, leads to a population of 100nm diameter unilamellar DOPA:OA vesicles.Again each vesicle is assumed spherical (Φ = 1) with aqueous volume Ω 0 , but now part of the bilayer consists of DOPA phospholipid in molecular fraction ρ.The molar concentration of DOPA:OA vesicles in the extruded suspension is approximately where the aggregation number N DOPA:OA is now calculated as the sum of both the OA lipids and DOPA phospholipids making up each closed bilayer.In turn, the OA lipid monomer concentration inside/outside the vesicles is [L] DOPA:OA eq , the CVC value for the model DOPA:OA vesicles, maintaining them at equilibrium.
Competition starts (t = 0) when the extruded vesicle solutions above are mixed.We mix a volume Ω p of each solution, creating a new mixed system of volume Ω e = 2Ω p , containing DOPA:OA vesicles in number N A Ω p C DOPA:OA ves and OA vesicles in number N A Ω p C OA ves .The initial lipid monomer concentration in the environment becomes 1 2 ([L] DOPA:OA eq + [L] OA eq ).Throughout mixing, and during competition, buffer concentration is constant at [B] in all solutions, at a value high enough for vesicles to maintain approximately constant volume Ω 0 .
Modelling the opposite scenario, namely a fixed population of pure OA vesicles mixed with a variable population of DOPA:OA vesicles, just requires switching the R multiplier from ( 17) to (18).

Mixing for Osmotically-Driven Competition
When a fixed population of isotonic OA vesicles is to be mixed in ratio R with a variable population of swelled OA vesicles, again two extruded vesicle suspensions are prepared.The first is prepared in buffer at molar concentration [B] and extruded through 100nm diameter pores, leading to unilamellar OA vesicles at Φ = 1 in molar concentration The second suspension is prepared in a solution which contains an additional membrane impermeable (or slowly permeating) solute, such as sucrose, mixed with the buffer, increasing the overall molar concentration of osmotically active species to [B] 0 ≥ [B] + 0.7.This suspension is made of unilamellar OA vesicles at Φ = 1 in molar concentration RC isotonic ves , and each vesicle encapsulates buffer at concentration [B] 0 .
The buffer concentration outside the vesicles in the second suspension is then reduced to [B], making the external solution hypotonic with respect to the vesicle interiors.The vesicles swell to maximum size, and then transiently break, allowing for the escape of buffer molecules in excess.They later reseal with a residual buffer gradient of [B] max ∆ = 0.16M across the membrane, corresponding to a maximum osmotic pressure of 4 atm (Chen et al., 2004).In our model, each vesicle is therefore assumed to swell to volume Ω = Ω 0 (1 + (0.16/[B])), which remains constant for the duration of competition.
The decrease of the environmental buffer concentration is considered to take place at the same instant of mixing with the initial isotonic population.This defines the initial condition (t = 0) when competition starts.The mixed overall volume is Ω e = 2Ω p where isotonic vesicles number N A Ω p C isotonic ves , and the swelled vesicles number an R multiple of this.The lipid monomer concentration in this new, larger environment is initially [L] OA eq .

Vesicle Breakage
During competition, to a first approximation, we assume that all of the original vesicles remain intact, with none breaking apart through excessive osmotic stress.By using the Morse equation for osmotic pressure and data supplied in (Chen et al., 2004, Supplementary Material), we were able to calculate an approximate burst tolerance ≈ 0.21 for our model pure oleate vesicles, where these vesicles burst through excessive osmotic pressure when Φ < 1 − .Pure OA vesicles reached a minimum of Φ = 0.77 in our phospholipid-driven competition simulations, and a minimum of Φ = 0.70 in our osmotic-driven competition simulations reported in Fig. 4.These values do not overly exceed the burst tolerance.

Control Experiments: Mixing With Buffer
Mixing a vesicle population with a buffer solution is modelled as doubling the current system volume and diluting the initial vesicle density to one half.In this case, we assume that the buffer solution contains no vesicles, but free lipid monomer at concentration just below the CVC of oleic acid2 .

Two Competing Populations: Comparison with Experimental Results
Figure 4 compares predictions made by our kinetic model against experimentally reported surface growth of vesicles in phospholipid-driven (Budin and Szostak, 2011) and osmotically-driven (Chen et al., 2004) competition.
Top figures 4A and 4B show the dynamics of surface area change in phospholipid-driven competition.Figure 4A details, in real time, the relative surface area of a tracked population of DOPA:OA (ρ = 0.1) vesicles, when this population is mixed 1 : 1 with (i) pure OA vesicles (green lines), (ii) DOPA:OA (ρ = 0.1) vesicles (blue lines) or (iii) buffer (black lines).In Fig. 4B, the tracked population is instead pure OA vesicles, which are mixed 1 : 1 with the same three options outlined above.
Stochastic simulation of our lipid kinetics model (Gillespie, 1976) correctly predicts that when mixed 1 : 1, DOPA:OA vesicles steal lipid and grow (green lines, 4A) at the expense of the pure OA vesicles, which shrink (blue lines, 4B).In this case, there is also fairly good quantitative agreement with the experimentally observed time courses, with the indirect effect alone (d = 0 lines) accounting for most of the surface area change in our model.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 phospholipid-driven 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.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 1.6 0.7 0.8 0.9 1 0.8 0.9 Gain surface Gain surface Gain surface

Lose surface
Lose surface Lose surface 5 Lipid competition tipping points.A Phospholipid competition between 30 model phospholipid-laden vesicles, each with DOPA fraction randomly assigned over the uniform interval 0 < ρ0 < 1. Depending on its initial DOPA fraction, each vesicle starts at a point on the horizontal blue line, and grows (green arrows) or shrinks (red arrows) to a point on the black line.The form of the black line is specific to this particular competing population, and is computed by ( 12).Orange crosses show agreement with equilibrated stochastic simulation of the model, validating the fast computation of competition equilibrium method.Competition 'tipping point' is shown by blue circle: any vesicle with ρ crit 0 > 0.573 gains lipids from its competitors.B Phospholipid competition in four different populations of 30 model vesicles, with DOPA fraction randomly assigned over uniform intervals (i) 0 < ρ0 < 0.25, (ii) 0 < ρ0 < 0.5, (iii) 0.25 < ρ0 < 0.75 and (iv) 0.3 < ρ0 < 1.0, demonstrating the context-dependence of competition.C Osmotic competition between 30 model oleate vesicles each swelled by extra internal sucrose, randomly assigned over the uniform interval 0 ≤ [B]∆ < 0.16 molar.Any vesicle starting at tension state Φ crit 0 > 0.802 gains lipids from its competitors.See text for full discussion.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.1) vesicles as more OA vesicles are added (Fig. 4C).In the opposite scenario, we also verified that the model shows a plateau in the shrinkage of pure OA vesicles as more DOPA:OA (ρ = 0.1) vesicles are added (Fig. 4D).In both cases, the indirect effect (d = 0 lines) alone drives the majority of the surface size change, with the direct effect then 'tuning' the fit to experimental outcomes.
Importantly, the general outcome of phospholipiddriven competition in our model is for vesicles stealing lipid to grow in surface and finish at high Φ > 1 values (excess surface, flaccid), and for vesicles losing lipid to suffer reduced surface, finishing at Φ < 1 values (osmotically tense, spherical).This is observed experimentally, and indeed provides the basis for the conjecture that phospholipid-laden vesicles are more likely to divide spontaneously when gentle external shearing forces are applied (Budin and Szostak, 2011, p5250).
Moving to osmotically-driven competition, Fig. 4E shows stochastic simulation of a swelled population of vesicles competing with an initially isotonic (non-swelled) population.Simulation outcomes match quite well the experimental best-fit time courses, in particular for the growth of the swelled vesicles (not so accurately for the shrinkage of the non-swelled vesicles).In any case, it must be noted that the original experimental data 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 non-swelled vesicles will cause the shrinkage of the non-swelled vesicles to plateau, rather than to continue (note the logarithmic scale of Fig. 4F).
As with phospholipid-driven competition, in our model, the pure OA vesicles involved in osmoticallydriven competition can finish at a variety of surface sizes.However, unlike phospholipid-driven competition, all OA vesicles will finish with the same Φ < 1 value (equal osmotic stress).Here, surface changes do not translate to final differences in Φ, partly because the vesicles start with different aqueous volumes.This residual osmotic swelling is also observed experimentally in vesicles stealing lipid through osmotically-driven competition.In fact, it stands as the main criticism of the osmotically-driven competition scenario: swelled vesicles have to overcome a stronger energetic barrier in order to divide, making this an improbable route to spontaneous vesicle division (Adamala and Szostak, 2013).
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 phospholipid-laden 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 population of phospholipid-laden 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 57%. 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 ( 12)), which includes the number and composition of competitor vesicles present3 .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 can drive growth (Budin and Szostak, 2011, p5251).
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 direct effect makes the competition tipping point slightly lower, but no general statement can be made about the extent to which it affects vesicle growth, for this again 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).
Figure 5C shows that in a heterogeneous population where pure OA model vesicles are swelled with resid-ual buffer up to 0.16M, vesicles with low initial Φ values steal lipid from those with higher (less swelled) Φ values, with the tipping point between growing and shrinking at Φ crit 0 = 0.8.As a last remark, orange crosses marked on Figs 5A and 5C show that full stochastic simulations of the model (run all the way to equilibrium) agree with and thus validate the fast computation of competition equilibrium method.

Theoretical Predictions Beyond Current Experimental Results
Finally, we were able to explore more widely some of the the parameter space for phospholipid-driven and osmotically-driven competition, using our model to make some predictions.Figure 6A shows the stoichiometry results of phospholipid-driven competition in this wider context.A population of DOPA:OA (ρ = 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 ρ pop2 0 = 0 shows the result reported in Fig. 4C as the red line.Interestingly, this figure predicts that the absolute growth of the second population of vesicles will be maximal, when their phospholipid fraction is around 50%, and will decline again towards no overall growth as the phospholipid fraction approaches 100%.Figure 6B explores the stoichiometry of osmotically-driven competition in a similar way to phospholipid-driven competition.A fixed population of swelled vesicles is mixed with a second population, and the trapped residual buffer inside vesicles in the second population, as well as the mixing ratio R to the second population, are varied.To conclude these predictions, Fig. 6C shows the effects of osmotically-driven versus phospholipid-driven 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.

IV. DISCUSSION
In this work, we have presented a theoretical model of the transfer kinetics of lipid molecules between vesicles.We have shown that reasonable rate equations chosen for simple lipid uptake and release allow to reproduce fairly well data from controlled laboratory experiments on phospholipid-driven competition and osmotically-driven competition.Furthermore, we have been able to predict the outcome of several yet-to-be-performed 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 phospholipid-driven competition is that the membrane  phospholipid fraction is approximately stationary with respect to the timescale of simple lipid transfer between the supramolecular structure (i.e., the membrane bilayer) and the aqueous solution (both inwards and outwards) 4 .In reality, the off-rate of a lipid molecule from a bilayer is inversely proportional to the number of carbon atoms in the acyl chain of the lipid concerned, and phospholipids do have a small non-zero transfer rate (with a half time from hours to days (Budin and Szostak, 2011)).If 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 number of membrane phospholipids P µ .With no remaining asymmetries in P µ to drive competition, all vesicles would finish with the same number of simple lipids L. The appearance and then disappearance of competition would follow the same type of dynamics as those experimentally reported for nervonic acid (Budin and Szostak, 2011, Fig. 4D, therein) which redistributes between vesicles (simulation results not shown).However, if vesicles contained a metabolism which synthesised phospholipid, then lasting P µ asymmetries between vesicles could be continually maintained as steady states, in spite of the exchanging P fraction.The results of this study can be interpreted as the competition advantage bestowed upon a vesicle by a membrane phospholipid fraction given that this fraction 4 Likewise, the assumption we make with osmotically-driven competition is that the residual buffer inside the vesicles permeates very slowly through the bilayer membrane.
is somehow maintained as constant.Not explicitly modelling phospholipid synthesis processes grants a simplified lipid scenario (i.e. a materially-closed system which subsequently settles to equilibrium) where some analysis can be carried out.
The next point that deserves discussion is the role of the direct effect in driving the phospholipid-driven competition simulations performed with our model.A first observation to make is that even when the direct effect is disabled (d = 0), the remaining indirect effect can account for the majority of the vesicle surface growth observed experimentally in phospholipid-driven competition (blue lines, Figs 4C and 4D).Thus, whilst a direct effect could improve the fit to experimental results, we should conclude from our treatment of the problem and the results obtained, that the indirect effect is the main mechanism driving vesicle growth dynamics.A second observation is that, as stated in the Results in reference to Fig. 5B, the exact contribution of the direct effect depends on the specific details of the competition scenario.In the case of the latter figure, the lipid release multiplier function r linearly decreases the simple lipid off-rate with increasing phospholipid content, but this 'context dependence' should also be true for different choices of function r.
One curiosity in the results (both in vitro and in silico) is how DOPA:OA (ρ = 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 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 [L] OA eq .As this happens, (20) implies that the DOPA:OA vesicles will be absorbing more and more L lipids, in order to grow to a size in equilibrium with the external lipid monomer concentration.The DOPA:OA growth is thus halted only by the number of lipids in the system being limited to L t .In our kinetics model, this continuous growth happens with or without the direct effect present.The direct effect can drive larger growths when vesicles are small, but as surface area increases and the direct effect diminishes, it is the indirect effect which persists and continues to drive growth as more OA vesicles are added.
A final point worth highlighting is that when the lipid uptake function u given in ( 5) 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.In this case, 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 [L] OA eq .Rather, the DOPA:OA have much slower growth, and they even have a finite stable size when the outside lipid monomer concentration is exactly [L] OA eq .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.
Understanding in full detail the dynamics of these colloidal systems is certainly not an easy task.In any case, we consider this work just as a step further in the development of semi-realistic, coarse-grained descriptions of phenomena that, in reality, are extremely complex.Self-assembly 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 (Mouritsen, 2005;Rasmussen et al., 2009), constitute a tremendous challenge, indeed, both for experimental and theoretical 'systems chemistry' research (Ruiz-Mirazo et al., 2014) and for synthetic biology (Solé, 2009;Solé et al., 2007).In particular, the connection between basic metabolic reaction networks and membrane dynamics (including stationary growth and division cycles (Mavelli and Ruiz-Mirazo, 2013)) 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.
FIG. 1Two mechanisms of phospholipid-driven growth.A Indirect effect, whereby the presence of phospholipid in a vesicle membrane drives growth simply through a geometric asymmetry: only the lipid section of the bilayer (grey) is able to release lipid (orange arrows) whereas the whole of the bilayer surface (made of lipids and phospholipids) is able to absorb lipid monomer (green arrows).Phospholipid fraction is pictured as one continuous block to highlight the principle only.The indirect effect can be created also by nonlipid surfactant molecules (e.g.peptides) residing long enough in the membrane to increase surface absorption area.B Direct effect, whereby the acyl tails of the phospholipids have high affinity for packing closer to each other and increasing bilayer order, thus making the exit of the simple lipids more difficult.The direct effect is specific to the molecular structure of phospholipids.In both cases, growth eventually stops when the phospholipid fraction in the membrane becomes diluted.

FIG. 4
FIG. 4 Comparison between kinetic model predictions and experimental results.Top plots show dynamics of phospholipid-driven competition.A Surface growth of DOPA:OA vesicles over time (green lines) and B surface shrinkage of OA vesicles over time (blue lines), when a population of DOPA:OA (ρ = 0.1) vesicles are mixed 1:1 with pure OA vesicles (following our mixing procedure with Ωp = Ω dyn , 25 DOPA:OA are mixed with 20 OA).Coloured dots in figure backgrounds reproduce original experimental results from (Budin and Szostak, 2011, Figs 1A, 1B therein) respectively.Middle plots show vesicle stoichiometry effects in phospholipid-driven competition.C Continued surface growth of DOPA:OA population as more OA vesicles added and D plateau in surface shrinkage of OA vesicles as more DOPA:OA vesicles added.Black markers in figure backgrounds reproduce experimental results from(Budin and Szostak, 2011, Figs 1C, 1D therein)  respectively.Bottom plots show osmotically-driven competition results.E Growth dynamics of swelled OA vesicles (blue line) and shrinkage of isotonic vesicles (red line) compared against experimental best-fit exponential decay curves (grey lines) from(Chen et al., 2004, Figs  1D, 1B therein)  respectively.F Stoichiometry effects in osmotically-driven competition: shrinkage of OA vesicle surface reaches a plateau as more swelled vesicles are added.Black markers in figure background reproduce experimental results from(Chen  et al., 2004, Fig. 2A therein).See text for discussion.

Figs
Figs 4C and 4D correspond to the surface sizes reached in the limit of time in Figs 4A and 4B respectively.

TABLE I
Vesicle competition model parameters.