# Wildfire Smoke Plume Segmentation Using
Geostationary Satellite Imagery

Jeff Wen  Marshall Burke

###### Abstract

Wildfires have increased in frequency and severity over the past two decades,
especially in the Western United States. Beyond physical infrastructure damage
caused by these wildfire events, researchers have increasingly identified
harmful impacts of particulate matter generated by wildfire smoke on
respiratory, cardiovascular, and cognitive health. This inference is difficult
due to the spatial and temporal uncertainty regarding how much particulate
matter is specifically attributable to wildfire smoke. One factor contributing
to this challenge is the reliance on manually drawn smoke plume annotations,
which are often noisy representations limited to the United States. This work
uses deep convolutional neural networks to segment smoke plumes from
geostationary satellite imagery. We compare the performance of predicted plume
segmentations versus the noisy annotations using causal inference methods to
estimate the amount of variation each explains in Environmental Protection
Agency (EPA) measured surface level particulate matter $<$2.5µ$m$ in diameter
($\textrm{PM}_{2.5}$).

remote sensing, segmentation, computer vision, wildfire smoke

##  1 Introduction

Since the early 1980s, scientists at the National Environmental Satellite,
Data and Information Service (NESDIS) have explored the use of environmental
satellites to detect fires while other researchers identified the potential of
using satellite imagery for wildfire smoke detection (McNamara et al., 2004;
Chung & Le, 1984; Svejkovsky, 1985). Although NESDIS implemented algorithmic
fire detection with the introduction of the Hazard Mapping System (HMS) in
2002, wildfire smoke is still manually annotated by analysts using multi-frame
animations of Geostationary Operational Environmental Satellite (GOES) imagery
(McNamara et al., 2004; Ruminski et al., 2007). We combine the methods being
developed in computer vision with geostationary satellite imagery to identify
wildfire smoke plumes in near real-time, providing a method to improve
analysis with more accurate identification of smoke and potentially extend
analysis beyond the US.

Researchers have applied machine learning techniques to related problems with
mixed results. Specifically, Wan et al. (2011) used unsupervised learning
approaches to cluster smoke in RGB color images using both a sharpening and
mixture model approach (Wan et al., 2011). Their exploratory results showed
that it was possible to identify smoke plumes by analyzing satellite imagery,
but the unsupervised approach limited opportunities to examine the out-of-
sample performance of the models. Wolters and Dean present two approaches to
smoke segmentation by combining labeled hyperspectral images with logistic
regression models (Wolters & Dean, 2015, 2017). Their more recent work
illustrates the ability of auto-logistic regression models to capture the
spatial association between neighboring pixels, which allows for smoothing
over predictions. However, while the results are positive, the model’s error
rates are still relatively high.

More recent literature illustrates the potential of using neural networks to
estimate local pollution levels (Li et al., 2017; Hong et al., 2019).
Similarly, researchers have designed custom neural network architectures for
classifying different classes of aerosol including smoke, but the
classification does not isolate and segment plumes, which makes it challenging
to quantify the contribution wildfire smoke versus other sources of aerosol in
downstream analysis (Ba et al., 2019). Filonenko et al. (2017) compares the
performance of different CNN architectures for smoke detection but does not
consider segmentation of plumes in the images (Filonenko et al., 2017).

Neural networks have been successfully applied in varied image segmentation
tasks such as street-level scene labeling, neural structure labeling in
biomedical electron microscopy images, road network segmentation, and aerial
image segmentation (Badrinarayanan et al., 2015; Ronneberger et al., 2015;
Mnih, 2013; Marmanis et al., 2016). These techniques have allowed applied
researchers to automate previously manual, time-intensive annotation tasks
with models that can be trained and adapted to different problems.

Ramasubramanian et al. (2019) use a convolutional neural network architecture
with 6 bands of GOES-16 imagery as input and HMS satellite-derived annotations
that were subsequently corrected by experts as labels (Ramasubramanian et al.,
2019). Their model uses a single timestamp as input and relies on subject
matter experts to manually correct the imperfect labels, which resulted in a
dataset size of approximately 120 scenes containing smoke. The authors note
that their model performs smoke pixel classification using the pixel of
interest in addition to an input neighborhood around the pixel of interest.
This approach results in predicted boundaries that extend beyond the visible
smoke boundary and relies on a manually defined input-pixel neighborhood size.
This along with the need for expert corrected labels makes it difficult to
scale the approach to further analysis.

Larsen et al. (2021) utilize a fully convolutional network (FCN) to identify
wildfire smoke plumes from satellite observations obtained by the Himawari-8
geostationary satellite situated over East Asia and Western Pacific regions
(Larsen et al., 2021). They use a deterministic cloud-masking algorithm to
generate smoke versus non-smoke pixel labels, which are then used to train the
FCN. While the approach is similar, the manually annotated smoke plumes in our
work may better capture variability in the visual representation of smoke
compared to a deterministic algorithm.

Other recent literature similarly use deep convolutional neural networks
(CNNs) for identifying smoke plumes from UAV and drone images as well as
synthetically generated smoke images (Frizzi et al., 2021; Yuan et al., 2019).
These works show promising performance but differ from the work presented here
as the aerial or synthetic images used provide a side-view angle of smoke
rather than a top-down view. Li et al. (2018) and Zhu et al. (2020) both use
3D convolution based CNN approaches to segment smoke plumes from videos, but
these video sequences also present a side-viewing angle (Li et al., 2018; Zhu
et al., 2020). While the side-view images may provide higher resolution, their
availability can be inconsistent compared to geostationary satellite
observations. Given the temporal frequency of geostationary satellite imagery,
future work can explore using video segmentation networks on sequences of
geostationary images.

Figure 1: First row: Model input, ground truth annotation label, and model
predicted segmentation after thresholding. Second row: Additional input
channel 07, channel 11, and MERRA-2.

In this work, we adapt CNNs to automatically identify wildfire smoke plumes
from satellite imagery. We use satellite observations as input and HMS
annotated smoke plumes as the target labels to train our models. Because
analysts often generate these plume annotations using multiple hours of
satellite imagery, we investigate different methods for improving training
with these noisy labels. To validate our models, we compare the predicted
smoke plume segmentations against the HMS annotated smoke plumes on PM2.5
measurements from EPA monitoring stations. Better understanding of the spatial
extent of wildfire smoke is especially important in downstream research, such
as in the environmental economics literature, where accurate identification of
smoke plumes could lead to better estimates of the causal effect of smoke
exposure.

##  2 Methodology

###  2.1 Data

The main imagery data used are satellite earth observations captured by the
Geostationary Operational Environmental Satellite (GOES-16) positioned over
the eastern United States. GOES-16 along with the west coast GOES-17 satellite
provide frequent (every 5 minutes for the Continental US) multi-spectral
measurements that are largely used for weather modeling and storm tracking
over the Atlantic and Pacific oceans respectively. Images are downloaded if
there are associated smoke plumes during 2018 from the HMS annotations. In
order to ensure that the imagery captures radiance from the earth’s surface,
images were limited to between 12PM and 11PM Coordinated Universal Time (UTC).
Furthermore, due to the large number of daily observations, we limited
training images to California and Nevada from the 2018 wildfire season between
May-December of 2018. However, we perform testing on unseen data between May-
October of 2019 and 2020. For testing on external data from 2019-2020, we
download 3 images per day between May-October. These multiple images allow us
to identify if multiple plumes are overhead on a given day. Given the
frequency of GOES imagery, in future work we hope to increase the temporal
resolution of the plumes.

Table 1: Input channels Name | Wavelength | Type | Primary uses
---|---|---|---
Blue |  0.47 $\mu$m | Visible | Aerosols
Red |  0.64 $\mu$m | Visible | Clouds & ash
Veggie |  0.86 $\mu$m | Near-IR | Vegetation
Shortwave |  3.9 $\mu$m | IR | Fire hotspots
Cloud-Top |  8.4 $\mu$m | IR | Cloud-top

After collecting the smoke plume data and the associated GOES-16 imagery from
the Amazon Web Services (AWS) NOAA open data repository, the raw imagery is
reprojected and a pseudo true-color composite is generated from the red, blue,
and ”green” bands using the SatPy package (Raspaud et al., 2021). Although the
GOES-16 satellite carries multiple sensors, it lacks a green band, which must
be generated to create a true-color composite (Bah et al., 2018).
Additionally, the 07 and 11 channels are also included in experiments as these
near-infrared channels can be used for fire hotspot detection and cloud-top
identification (examples in Figure 1). Each of the downloaded images was saved
as a 1200x1200 .geotiff image for a total of 615 images. These were then
randomly cropped to generate up to 15 different 300x300 images where 60% of
the crops had a smoke plume to deal with class imbalance. As a result, there
were 6825 images which was split into 70% training, 15% validation, and 15%
testing. We further ensured that the cropped images generated from the same
base image could only be in one set of data to reduce potential data leakage.

We pair the satellite imagery with smoke plume labels from the HMS archives.
While the labels extend back to 2006, we use observations from the most recent
generation of geostationary satellite (GOES-16) launched in mid-2017 as this
subset of data allows for more spectral bands to be used in the analysis.
Analysts often use multi-hour animations to draw the extent of the smoke
plumes, which presents challenges as the extent over time might not match the
boundaries of smoke plumes in a single image due to wind or other
meteorological factors (Brey et al., 2018). Although the current labels are
manually annotated from satellite imagery, one goal of this research is to
identify the extent to which we can learn meaningful segmentations from these
imprecise labels. We attempt to improve training on noisy labels by providing
additional signal using the Modern-Era Retrospective analysis for Research and
Applications, Version 2 (MERRA-2) aersol optical thickness (AOT) measurements
of particulate matter as an additional input channel (Gelaro et al., 2017).
The MERRA-2 data combines multiple sources of aerosol optical density
information mainly derived from the Moderate Resolution Imaging
Spectroradiometer (MODIS) instrument aboard NASA’s Terra and Aqua satellites.
We discuss 2 additional methods to account for the noisy labels in the next
section.

In the rest of the paper, we will refer to the true-color image as 1 band and
the true-color, channel 07, and channel 11 combined as 3 bands. Lastly, 4
bands refers to using the true-color, channel 07, channel 11, and the MERRA-2
in the input.

###  2.2 Model

Our baseline experiments utilize an adapted U-Net neural network architecture,
originally designed for biomedical image segmentation, to segment smoke plumes
from a pseudo true-color RGB image (Ronneberger et al., 2015). The adapted
network keeps the same number of layers in the encoding and upsampling blocks
but decreases the number of convolutional filters. Furthermore, we replace the
ReLU activation function with PReLU activations for improved stability in
training. We also include multiple spectral bands in subsequent experiments as
these bands (shown in Table 1) provide contextual information of fire hotspot
location and cloud-top phase in channel 07 and 11 respectively.

Across the model training, we use the Adam optimization with a learning rate
starting at $\alpha$=5e-5 and stepping down by $\gamma$=0.1 every 9 epochs
(Kingma & Ba, 2017). Each model was trained for a total of 21 epochs using a
batch size of 16. We compared binary cross entropy (BCE) and mean absolute
error (MAE) losses during the training and validation process. To track
accuracy during training and validation, we kept track of the average Dice
coefficient.

$\text{Dice Coefficient}=\frac{2|A\cap B|}{|A|+|B|}$ (1)

This statistical similarity metric first proposed by Lee Raymond Dice in 1945
and shown in Equation 1 calculates 2 times the amount of overlapping pixels
between the predicted ($A$) and ground truth mask ($B$) divided by the total
number of pixels in both masks (Dice, 1945).

Additionally, we consider two methods for explicitly handling noisy labels in
the training process. Specifically, we experiment with using mean absolute
value of error (MAE) loss, which has been shown to be tolerant to label noise,
and data sampling to specify that only low loss training samples contribute to
the gradient updates (Ghosh et al., 2017; Xue et al., 2019). The loss sampling
strategy allows the model to run the prediction on a batch of inputs as usual,
but before calculating the average loss across samples, we identify the
training sample in the batch that produced the highest loss. Then, we set the
loss for that example to 0 so that it does not contribute to making weight
updates. As mentioned by Xue et al., this approach assumes that as the model
performance improves, particularly noisy samples can result in high loss,
which can have large impact on the weight updates. This training strategy
mitigates the effect of these samples.

Figure 2: The model predicted segmentation on the right more closely matches
the visible smoke in the input imagery compared to the HMS annotated smoke
plumes in the middle, which cover most of California. The input on the left is
a test image from September 8th, 2019 that was not used for training or
validation.

###  2.3 External validation

While the DICE coefficient is used for validation during training, we apply
the models on unseen images from 2019 and 2020 to estimate the performance of
these models compared to the hand annotated HMS smoke plumes from the same
years. This allows us to test the models in a setting most similar to
downstream inference tasks where the models would be used to identify smoke
plumes across a certain time frame, then used for additional analysis.

Specifically, we leverage econometric tools for causally identifying the
contribution of wildfire smoke to changes in ground level PM2.5 readings as
measured by the EPA. We use a quasi-experimental approach to exploit variation
in smoke and PM2.5 over time to estimate the effect of identified wildfire
smoke plumes on PM2.5 readings. We use fixed effect panel regressions (shown
in Equation 2) with station fixed effects to account for time-invariant
unobserved effects such as the fact that different stations may have
unobserved characteristics correlated with smoke exposure. We set the station
PM2.5 reading as the dependent variable and consider the HMS annotations as
well as different model predicted smoke plumes as the independent variable.

$\text{PM2.5}_{it}=\beta_{1}\text{Smoke}_{it}+\alpha_{i}+\epsilon_{it}$ (2)

This approach allows us to compare performance even in the presence of noisy
annotation labels by measuring performance against an external ground truth.
In the above equation, the $\alpha_{i}$ captures the EPA station fixed
effects, which would handle time-invariant station level differences. The
$\text{Smoke}_{it}$ variable is determined by the presence or absence of smoke
plumes overlapping an EPA station where the smoke plumes come from either the
HMS annotations or model predicted segmentations.

Table 2: Model performance on validation set. Bold rows represent the models that achieve the best validation DICE coefficient for a specific loss function. The asterisk denotes that the model weights were updated after removing the highest loss sample per batch. Bands | Loss | Avg. Loss | DICE
---|---|---|---
1 | BCE | 0.2535 | 0.0948
3 | BCE | 0.2236 | 0.1074
3* | BCE | 0.2313 | 0.1008
4 | BCE | 0.1884 | 0.1028
1 | MAE | 0.0986 | 0.2635
3 | MAE | 0.0986 | 0.2649
3* | MAE | 0.0986 | 0.2655
4 | MAE | 0.0986 | 0.2655

To compare the HMS annotated smoke plumes against the model predicted
segmentations, we consider the adjusted R2 (Adj.R2) and within-adjusted R2 (W
Adj.R2) metrics. The Adj.R2 metric measures the total amount of variation in
PM2.5 measured at EPA monitoring stations explained by the smoke plumes and
the W Adj.R2 metric represents the amount of variation in PM2.5 that is
explained ”within” the PM2.5 station unit by the different sources of smoke
data. This distinction is important and we prefer the ”within” metric because
it is a better indicator of ability to explain variation in PM2.5 using the
different smoke data as it only considers remaining variation separate from
time-invariant unobserved station level differences.

##  3 Results

Our preliminary results show that qualitatively the BCE loss models are able
to learn to segment smoke from the input imagery. The top right panel in
Figure 1 shows the performance of the model on a validation sample and the
rightmost panel in Figure 2 and Figure 3 show segmentations of the model
trained with BCE loss and 3 input channels on previously unseen test images.
The segmentations appear to more precisely locate the smoke in the input
imagery on the left compared to the annotations in the middle. Even though
there is cloud cover in Figure 1, the model is able to learn to differentiate
between the cloud cover and the smoke plumes. Further quantitative results are
shown in Table 2, where the bolded rows represent the best validation DICE
coefficient when using either of the loss functions during training. The
asterisk denotes that the model weights were updated only after removing
training samples that had the highest loss per batch (Xue et al., 2019).

Although the MAE loss resulted in higher overall accuracy, the qualitative
results indicated that these models were predicting ”no-smoke” for nearly all
pixels and were unable to deal with the label imbalance in the images. Even
with the addition of the additional step to sample only low loss training
samples, the model performance with MAE loss was still qualitatively unable to
identify smoke pixels compared to the models trained with BCE loss.

Table 3: HMS annotated vs. model predicted plume performance on EPA PM2.5 measurements. Different sources of smoke plumes are validated against external PM2.5 measurements using a fixed effect panel regression approach to compare variation explained. The asterisk denotes the loss sampling model results. | Adj. R2 | W Adj. R2
---|---|---
Annotations | 0.2316 | 0.1946
1 band | 0.2707 | 0.2367
3 band | 0.2715 | 0.2374
4 band | 0.3038 | 0.2703
3 band* | 0.2987 | 0.2653

When we evaluate the model performance on previously unseen data and compare
the performance against the annotations, we see that the model predictions are
able to explain more of the variation in surface level PM2.5 readings
regardless of model specification (results shown in Table 3). The higher
values for the model predicted smoke data suggest that while these results are
preliminary, the CNN segmented smoke plumes may explain more of the changes in
PM2.5 compared to the hand annotated smoke plumes. It is important to note
that the hand annotations are meant to capture visible smoke, while the
particulate matter captured at EPA stations might not be visible from
satellite imagery. Therefore, the models leveraging the MERRA-2 channel as
additional input might not be a fair comparison against the hand annotations,
which are mainly generated using visible imagery. However, as seen in Table 3,
even the model trained only using RGB as input appears to capture more of the
”within” station variation than the hand annotations.

While in Table 2 the model trained by removing the influence of the highest
loss sample per batch appeared to perform worse than the best 3-band model
using BCE loss, in the external validation (shown in Table 3), this model was
able to explain more of the total and within variation in PM2.5. Furthermore,
this model was nearly able to match the performance of the 4-band model
illustrating that explicitly ignoring high loss training samples when there
are noisy labels could be beneficial for downstream tasks. This might not be
apparent from the validation accuracy metrics because the validation dataset
also includes noisy labels, which make it challenging to gauge model
performance. This underscores the importance of clean validation data
especially in the presence of noisy labels.

Figure 3: Additional examples show model predictions on unseen test images
with a lot of smoke in the first row and no smoke in the second row. Model
input, ground truth annotation label, and model predictions are displayed in
both rows from left to right.

##  4 Conclusion

As wildfires continue to worsen, it is increasingly important to quantify the
effects of wildfire smoke exposure on society. In this work, we used an
adapted U-Net architecture to segment smoke plumes from geostationary
satellite imagery with the goal of improving our understanding of the spatial
extent of smoke. We further leveraged quasi-experimental methods to compare
the variation in EPA station PM2.5 measurements that could be explained by
either model predicted segmentations or hand annotated plumes. While smoke
plumes have been manually annotated in the United States since the 2000s, our
results suggest that automated segmentation methods are at least qualitatively
comparable to the annotated smoke plumes and explain more of the within
station PM2.5 variation in external validation. These early results show the
potential of adapting neural networks for improving downstream inference even
with noisy labels.

In future work, we look to focus on additional methods to improve the model
robustness to noisy labels. We also hope to extend and develop these methods
to identify smoke across the globe and over time to better inform the impacts
of wildfires.

### Software and Data

The Geostationary Operational Environmental Satellite (GOES) data used in this
work is available for download at https://registry.opendata.aws/noaa-goes.
Hazard Mapping System smoke data is available at
https://www.ospo.noaa.gov/Products/land. The models were developed using
PyTorch version 1.7.1 and the fixed effects estimators were estimated using
the Fixest R package version 0.8.4 (Paszke et al., 2019; Bergé, 2018). The
models and data will be made available at the time of publication.

### Acknowledgements

We thank reviewers whose comments and suggestions helped clarify and improve
this manuscript. We also thank the Sustainability and Artificial Intelligence
Lab for discussions.

## References

  * Ba et al. (2019) Ba, R., Chen, C., Yuan, J., Song, W., and Lo, S.  Smokenet: Satellite smoke scene detection using convolutional neural network with spatial and channel-wise attention.  _Remote Sensing_ , 11:1702, 07 2019.  doi: 10.3390/rs11141702.
  * Badrinarayanan et al. (2015) Badrinarayanan, V., Kendall, A., and Cipolla, R.  Segnet: A deep convolutional encoder-decoder architecture for image segmentation.  _CoRR_ , abs/1511.00561, 2015.  URL http://arxiv.org/abs/1511.00561.
  * Bah et al. (2018) Bah, M., Gunshor, M., and Schmit, T.  Generation of goes-16 true color imagery without a green band.  _Earth and Space Science_ , 5(9):549–558, 2018\.
  * Bergé (2018) Bergé, L.  Efficient estimation of maximum likelihood models with multiple fixed-effects: the R package FENmlm.  _CREA Discussion Papers_ , (13), 2018.
  * Brey et al. (2018) Brey, S. J., Ruminski, M., Atwood, S. A., and Fischer, E. V.  Connecting smoke plumes to sources using hazard mapping system (hms) smoke and fire location data over north america.  _Atmospheric Chemistry and Physics_ , 18(3):1745–1761, 2018.  doi: 10.5194/acp-18-1745-2018.  URL https://www.atmos-chem-phys.net/18/1745/2018/.
  * Chung & Le (1984) Chung, Y.-S. and Le, H.  Detection of forest-fire smoke plumes by satellite imagery.  _Atmospheric Environment (1967)_ , 18(10):2143–2151, 1984.
  * Dice (1945) Dice, L. R.  Measures of the amount of ecologic association between species.  _Ecology_ , 26(3):297–302, 1945.
  * Filonenko et al. (2017) Filonenko, A., Kurnianggoro, L., and Jo, K.-H.  Comparative study of modern convolutional neural networks for smoke detection on image data.  In _2017 10th International Conference on Human System Interactions (HSI)_ , pp. 64–68. IEEE, 2017.
  * Frizzi et al. (2021) Frizzi, S., Bouchouicha, M., Ginoux, J.-M., Moreau, E., and Sayadi, M.  Convolutional neural network for smoke and fire semantic segmentation.  _IET Image Processing_ , 15(3):634–647, 2021\.
  * Gelaro et al. (2017) Gelaro, R., McCarty, W., Suarez, M. J., Todling, R., Molod, A., Takacs, L., Randles, C. A., Darmenov, A., Bosilovich, M. G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., da Silva, A. M., Gu, W., Kim, G.-K., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J. E., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert, S. D., Sienkiewicz, M., and Zhao, B.  The modern-era retrospective analysis for research and applications, version 2 (merra-2).  _Journal of Climate_ , 30(14):5419 – 5454, 2017\.  doi: 10.1175/JCLI-D-16-0758.1.  URL https://journals.ametsoc.org/view/journals/clim/30/14/jcli-d-16-0758.1.xml.
  * Ghosh et al. (2017) Ghosh, A., Kumar, H., and Sastry, P. S.  Robust loss functions under label noise for deep neural networks, 2017\.
  * Hong et al. (2019) Hong, K. Y., Pinheiro, P. O., and Weichenthal, S.  Predicting global variations in outdoor pm2.5 concentrations using satellite images and deep convolutional neural networks, 2019.
  * Kingma & Ba (2017) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization, 2017.
  * Larsen et al. (2021) Larsen, A., Hanigan, I., Reich, B. J., Qin, Y., Cope, M., Morgan, G., and Rappold, A. G.  A deep learning approach to identify smoke plumes in satellite imagery in near-real time for health risk communication.  _Journal of exposure science & environmental epidemiology_, 31(1):170–176, 2021.
  * Li et al. (2017) Li, T., Shen, H., Yuan, Q., Zhang, X., and Zhang, L.  Estimating ground-level pm2.5 by fusing satellite and station observations: A geo-intelligent deep learning approach.  _Geophysical Research Letters_ , 44(23):11,985–11,993, Dec 2017.  ISSN 0094-8276.  doi: 10.1002/2017gl075710.  URL http://dx.doi.org/10.1002/2017GL075710.
  * Li et al. (2018) Li, X., Chen, Z., Wu, Q. J., and Liu, C.  3d parallel fully convolutional networks for real-time video wildfire smoke detection.  _IEEE Transactions on Circuits and Systems for Video Technology_ , 30(1):89–103, 2018.
  * Marmanis et al. (2016) Marmanis, D., Wegner, J. D., Galliani, S., Schindler, K., Datcu, M., and Stilla, U.  Semantic segmentation of aerial images with an ensemble of cnss.  _ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2016_ , 3:473–480, 2016.
  * McNamara et al. (2004) McNamara, D., Stephens, G., Ruminski, M., and Kasheta, T.  The hazard mapping system (hms)–noaa multi-sensor fire and smoke detection program using environmental satellites.  In _13th Conf. Satellite Meteorology and Oceanography_ , 2004.
  * Mnih (2013) Mnih, V.  _Machine Learning for Aerial Image Labeling_.  PhD thesis, University of Toronto, 2013.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S.  Pytorch: An imperative style, high-performance deep learning library.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. Curran Associates, Inc., 2019.  URL http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * Ramasubramanian et al. (2019) Ramasubramanian, M., Kaulfus, A., Maskey, M., Ramachandran, R., Gurung, I., Freitag, B., and Christopher, S.  Pixel level smoke detection model with deep neural network.  In _Image and Signal Processing for Remote Sensing XXV_ , volume 11155, pp. 1115515. International Society for Optics and Photonics, 2019.
  * Raspaud et al. (2021) Raspaud, M., Hoese, D., Lahtinen, P., Finkensieper, S., Holl, G., Dybbroe, A., Proud, S., Meraner, A., Zhang, X., Joro, S., joleenf, Roberts, W., Ørum Rasmussen, L., Méndez, J. H. B., Zhu, Y., Daruwala, R., strandgren, BENR0, Jasmin, T., Barnie, T., Sigurosson, E., R.K.Garcia, Leppelt, T., ColinDuff, Egede, U., LTMeyer, Itkin, M., Goodson, R., jkotro, and peters77.  pytroll/satpy: Version 0.27.0 (2021/03/26), March 2021.  URL https://doi.org/10.5281/zenodo.4638572.
  * Ronneberger et al. (2015) Ronneberger, O., Fischer, P., and Brox, T.  U-net: Convolutional networks for biomedical image segmentation.  _CoRR_ , abs/1505.04597, 2015.  URL http://arxiv.org/abs/1505.04597.
  * Ruminski et al. (2007) Ruminski, M., Kondragunta, S., Draxler, R., and Rolph, G.  Use of environmental satellite imagery for smoke depiction and transport model initialization.  In _16th Annual International Emission Inventory Conf.: Emission Inventories—Integration, Analysis, and Communications_ , 2007.
  * Svejkovsky (1985) Svejkovsky, J.  Santa ana airflow observed from wildfire smoke patterns in satellite imagery.  _Monthly Weather Review_ , 113(5):902–906, 1985\.
  * Wan et al. (2011) Wan, V., Braun, W. J., Dean, C., and Henderson, S. B.  A comparison of classification algorithms for the identification of smoke plumes from satellite images.  _Statistical Methods in Medical Research_ , 20(2):131–156, 2011.  doi: 10.1177/0962280210372454.  URL https://doi.org/10.1177/0962280210372454.  PMID: 20889573.
  * Wolters & Dean (2015) Wolters, M. A. and Dean, C.  Issues in the identification of smoke in hyperspectral satellite imagery — a machine learning approach.  2015\.
  * Wolters & Dean (2017) Wolters, M. A. and Dean, C.  Classification of large-scale remote sensing images for automatic identification of health hazards.  _Statistics in Biosciences_ , 9(2):622–645, Dec 2017.  ISSN 1867-1772.  doi: 10.1007/s12561-016-9185-5.  URL https://doi.org/10.1007/s12561-016-9185-5.
  * Xue et al. (2019) Xue, C., Dou, Q., Shi, X., Chen, H., and Heng, P. A.  Robust learning at noisy labeled medical images: Applied to skin lesion classification, 2019.
  * Yuan et al. (2019) Yuan, F., Zhang, L., Xia, X., Wan, B., Huang, Q., and Li, X.  Deep smoke segmentation.  _Neurocomputing_ , 357:248–260, 2019.
  * Zhu et al. (2020) Zhu, G., Chen, Z., Liu, C., Rong, X., and He, W.  3d video semantic segmentation for wildfire smoke.  _Machine Vision and Applications_ , 31(6):1–10, 2020.


# Modelling Lipid Competition Dynamics in Heterogeneous Protocell Populations

Ben Shirt-Ediss  ICREA-Complex Systems Lab, Institut de Biologia Evolutiva,
CSIC-UPF, Barcelona, Spain  Logic and Philosophy of Science Department,
University of The Basque Country, Spain  Kepa Ruiz-Mirazo  Biophysics Unit
(CSIC-UPV/EHU), University of The Basque Country, Spain  Logic and Philosophy
of Science Department, University of The Basque Country, Spain  Fabio Mavelli
Chemistry Department, University of Bari, Italy  Ricard V.
Solé111Corresponding author<EMAIL_ADDRESS>ICREA-Complex Systems Lab,
Institut de Biologia Evolutiva, CSIC-UPF, Barcelona, Spain  Santa Fe
Institute, 1399 Hyde Park Road, Santa Fe NM 87501, USA

###### Abstract

In addressing the origins of Darwinian evolution, recent experimental work has
been focussed on the discovery of simple physical effects which would provide
a relevant selective advantage to protocells competing with each other for a
limited supply of lipid. In particular, data coming from Szostak’s lab suggest
that the transition from simple prebiotically plausible lipid membranes to
more complex and heterogeneous ones, closer to real biomembranes, may have
been driven by changes in the fluidity of the membrane and its affinity for
the available amphiphilic compound, which in turn would involve changes in
vesicle growth dynamics. Earlier work from the same group reported
osmotically-driven competition effects, whereby swelled vesicles grow at the
expense of isotonic ones. In this paper, we try to expand on these
experimental studies by providing a simple mathematical model of a population
of competing vesicles, studied at the level of lipid kinetics. In silico
simulations of the model are able to reproduce qualitatively and often
quantitatively the experimentally reported competition effects in both
scenarios. We also develop a method for numerically solving the equilibrium of
a population of competing model vesicles, which is quite general and
applicable to different vesicle kinetics schemes.

Protocells, chemical evolution, phospholipids, selection, self-assembly,
growth kinetics

##  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 (Smith and Szathmáry,, 1995; Luisi,,
2006). 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 ‘less-organised’ 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, aggregate-level 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 (Varela et al.,, 1974; Gánti,, 1975;
Dyson,, 1985; Segré et al.,, 2001; Solé et al.,, 2007; Macía and Solé,, 2007;
Mavelli and Ruiz-Mirazo,, 2007; Ruiz-Mirazo et al.,, 2011). 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 (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 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.

Figure 1:  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 non-lipid 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.

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

##  II Methods

###  II.1 Theoretical Model of Vesicle Competition

The competition model involves a set of $n$ vesicles

${\cal V}=\\{{\cal V}_{1},...,{\cal V}_{n}\\}$

each one characterized by a quadruple

${\cal V}_{j}=(L_{\mu}^{j},P_{\mu}^{j},L_{c}^{j},B_{c}^{j})$

embedded in a finite volume environment $\cal E$ defined by a triple
$(\Omega_{e},L_{e},B_{e})$. Each competing vesicle consists of a unilamellar
membrane of two different lipid types: a fixed number of phospholipids
$P_{\mu}$ (e.g. di-oleoyl-phosphatidic acid, DOPA) and a variable number of
simple fatty acid lipids $L_{\mu}$ (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 $\cal 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 $\cal E$ with constant
number $B_{e}$.

Vesicles compete with each other by consequence of uptaking/releasing simple
lipids $L$ from/to $\cal 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).

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

More precisely, each vesicle ${\cal V}_{j}$ is considered to release lipids to
both aqueous phases (at each side of the bilayer) at the equal rate of
$\lambda_{j}^{\text{out}}=k_{\text{out}}L_{\mu}^{j}{\textbf{r}}(\rho_{j})$,
and absorb lipids from each phase at rate
$\lambda_{j}^{\text{in}}=k_{\text{in}}S_{\mu}^{j}[L]{\textbf{u}}(\Phi_{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]_{c}^{j}=[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:

$L_{e}+\sum_{j=1}^{n}\left(L_{c}^{j}+L_{\mu}^{j}\right)-L_{t}=0$ (1)

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:

$\frac{dL_{c}^{j}}{dt}=k_{\text{out}}L_{\mu}^{j}{\textbf{r}}(\rho_{j})-k_{\text{in}}S_{\mu}^{j}[L]_{c}^{j}{\textbf{u}}(\Phi_{j})$
(2)

$\frac{dL_{\mu}^{j}}{dt}=-2k_{\text{out}}L_{\mu}^{j}{\textbf{r}}(\rho_{j})+k_{\text{in}}S_{\mu}^{j}([L]_{c}^{j}+[L]_{e}){\textbf{u}}(\Phi_{j})$
(3)

and $L_{e}$ can be deduced from constraint (1), once all $L_{c}$ and $L_{\mu}$
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_{\text{out}}$ of disassociating from the membrane, and function r modifies
this probability based on the current molecular fraction of phospholipid in
the membrane $\rho=\frac{P_{\mu}}{P_{\mu}+L_{\mu}}$. In order to account for
the _direct effect_ , we define function $0\leq\textbf{r}(\rho)\leq 1$ to be
monotonically decreasing with increasing $\rho$, 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:

$\textbf{r}(\rho)=1-d\rho$ (4)

where parameter $0\leq d\leq 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_{\mu}$ and function u, based on the dimensionless
geometric factor $\Phi=S_{\mu}/\sqrt[3]{36\pi\Omega^{2}}$ (where $\Omega$ is
vesicle aqueous volume, in the same units as $S_{\mu}$). We define a
conditional function

$\textbf{u}({\Phi})=\begin{cases}exp\left(\frac{1}{\Phi}-1\right),&\Phi<1\\\
1,&\Phi\geq 1\end{cases}$ (5)

to denote that lipid uptake is only increased when the the bilayer is stressed
($\Phi<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_{\mu}=\frac{1}{2}(L_{\mu}\alpha_{L}+P_{\mu}\alpha_{P})$, is the water-
exposed area of the inner bilayer leaflet, or alternatively the water-exposed
area of the outer bilayer leaflet of the vesicle. Membrane thickness is
therefore considered negligible. Uptake and release kinetic constants
$k_{\text{in}}$ and $k_{\text{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]_{\text{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, $\cal E$ has a fixed volume of
$\Omega_{e}$ litres, and each vesicle ${\cal V}_{j}$ has, in principle, a
variable volume internal water pool of
$\Omega_{j}=\Omega_{e}(L_{c}^{j}+B_{c}^{j})/(L_{e}+B_{e})$ litres. This
condition ensures that, at all times, the interior of each vesicle is isotonic
with respect to $\cal 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
$\Omega_{j}\approx\Omega_{e}(B_{c}^{j}/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 1.

###  II.2 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 associated to the full model are ignored
in order to keep only the logic of the problem (Fig. 3).

Figure 3:  Meanfield model of vesicle population dynamics. Considering
vesicles as simplified aggregates permits some analytical treatment.

The first simplification will be to ignore the internal structure of the
vesicles, describing them instead as coarse-grained ‘aggregates’, denoted by
pairs ${\cal 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
(${\textbf{u}}(\Phi)=1$) and independent of membrane phospholipid fraction
(${\textbf{r}}(\rho)=1$) respectively. Thus, the ODE system reduces to $n$
simplified equations, where for each aggregate:

$\frac{dL_{j}}{dt}=-k_{\text{out}}L_{j}+\frac{1}{2}k_{\text{in}}(L_{j}\alpha_{L}+P_{j}\alpha_{P})[L]_{e}$
(6)

Under these conditions, at equilibrium, the molar lipid concentration in the
environment $[L]_{e}=[L]_{\text{eq}}$ is related to the number of lipids and
phospholipids in an aggregate by the following function:

$f(L_{j},P_{j})=[L]_{\text{eq}}=\frac{2k_{\text{out}}}{k_{\text{in}}}\frac{L_{j}}{L_{j}\alpha_{L}+P_{j}\alpha_{P}}$
(7)

For a fixed number of phospholipids $P_{j}>0$, the mapping
$f:L_{j}\rightarrow[L]_{\text{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:

$\frac{dL_{i}}{dt}=-k_{\text{out}}L_{i}+\eta(L_{i}\alpha_{L}+P_{i}\alpha_{P})(L_{t}-\sum_{m=1}^{n}L_{m})$

$\frac{dL_{j}}{dt}=-k_{\text{out}}L_{j}+\eta(L_{j}\alpha_{L}+P_{j}\alpha_{P})(L_{t}-\sum_{m=1}^{n}L_{m})$

where $\eta={k_{\text{in}}}/{2N_{A}\Omega_{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

$\frac{L_{i}\alpha_{L}+P_{i}\alpha_{P}}{L_{j}\alpha_{L}+P_{j}\alpha_{P}}=\frac{L_{i}}{L_{j}}$
(8)

which leads to the following proportionality relation at equilibrium:

$L_{i}=\left(\frac{P_{i}}{P_{j}}\right)L_{j}$ (9)

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}\neq 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.

###  II.3 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_{\mu}\rightarrow[L]_{\text{eq}}$, like
(7), which gives the inside/outside lipid monomer concentration
$[L]_{\text{eq}}$ necessary to maintain a particular vesicle ${\cal V}_{j}$ at
equilibrium, given that this vesicle has a specific number of
lipids/phospholipids in the membrane, and a specific volume:

$f(L_{\mu}^{j},P_{\mu}^{j},\Omega_{j})=[L]_{\text{eq}}^{j}=\frac{2k_{\text{out}}}{k_{\text{in}}}\frac{L_{\mu}^{j}}{L_{\mu}^{j}\alpha_{L}+P_{\mu}^{j}\alpha_{P}}\frac{\textbf{r}(\rho_{j})}{\textbf{u}({\Phi_{j}})}$
(10)

The inverse of this function yields useful information: it is the mapping of
$[L]_{\text{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_{\mu}$ 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 ${\cal V}_{j}$ by using $f$
and binary searching for an $L_{\mu}$ which satisfies:

$f^{-1}([L]_{\text{eq}},P_{\mu}^{j},\Omega_{j})=L_{\mu}^{j}\;|\;f(L_{\mu}^{j},P_{\mu}^{j},\Omega_{j})-[L]_{\text{eq}}=0$
(11)

using appropriate search bounds (normally: $L_{\mu}^{\text{min}}=0$,
$L_{\mu}^{\text{max}}=L_{t}$).

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]_{\text{eq}}$. For each $[L]_{\text{eq}}$, we know that
each vesicle has a _unique_ number of membrane lipids $L_{\mu}$, because
$f^{-1}$ is itself one-to-one. This means that a certain $[L]_{\text{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):

$[L]^{*}=[L]_{\text{eq}}\;|\;\sum_{j=1}^{n}f^{-1}([L]_{\text{eq}},P_{\mu}^{j},\Omega_{j})+\\\
[L]_{\text{eq}}N_{A}\Omega_{e}-L_{t}=0$ (12)

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]_{\text{eq}}$ between
appropriate bounds, normally
$[L]_{\text{eq}}^{\text{(min)}}=\text{max}(f(L_{\mu}=0))$ over all vesicles,
$[L]_{\text{eq}}\text{(max)}=\text{min}(f(L_{\mu}=L_{t}))$ over all vesicles.
Finally, knowing $[L]^{*}$ allows to fully reconstruct the final sizes of all
vesicles at equilibrium by substituting $[L]_{\text{eq}}=[L]^{*}$ into (11)
for each vesicle ${\cal 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_{\mu}^{\text{crit}}$ separating those vesicles which have
lost lipid from those which have gained lipid, and is found by:

$P_{\mu}^{\text{crit}}=P_{\mu}\;|\;f^{-1}([L]^{*},P_{\mu}^{j},\Omega_{j})-\frac{2S_{\mu}-P_{\mu}\alpha_{P}}{\alpha_{L}}=0$
(13)

where $S_{\mu}=S_{\mu}^{0}$, again solvable by binary searching, this time in
the range $0\leq P_{\mu}\leq 2S_{\mu}^{0}/\alpha_{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_{\mu}$ to the value
desired.

The critical phospholipid number can be stated more usefully as the critical
phospholipid molecular fraction

$\rho_{0}^{\text{crit}}=\frac{P_{\mu}^{\text{crit}}\alpha_{L}}{2S_{\mu}^{0}+P_{\mu}^{\text{crit}}(\alpha_{L}-\alpha_{P})}$
(14)

a vesicle has in the initial condition, a time when all vesicles have a
surface of $S_{\mu}^{0}$. For osmotically-driven competition, the critical
volume separating shrinking vesicles from growing vesicles is found by
searching (13) for vesicle volume instead:

$\Omega^{\text{crit}}=\Omega\;|\;f^{-1}([L]^{*},P_{\mu},\Omega)-\frac{2S_{\mu}-P_{\mu}\alpha_{P}}{\alpha_{L}}=0$
(15)

where $S_{\mu}=S_{\mu}^{0}$. This may be alternatively stated as the critical
$\Phi$ in the initial condition:

$\Phi_{0}^{\text{crit}}=\frac{S_{\mu}^{0}}{\sqrt[3]{36\pi(\Omega^{\text{crit}})^{2}}}$
(16)

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 answer222Model 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 more
arguments to function $f$: the first would detail the head area for the
phospholipid in vesicle ${\cal 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)..

###  II.4 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 $\Omega_{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 $\Omega_{p}=\Omega_{\text{stoi}}$ litres (Table 1) 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 $\Omega_{p}=\Omega_{\text{dyn}}$ litres to be used
(to give $\tau$ 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 $\rho$ 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}\gg\text{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. Each vesicle is
assumed spherical ($\Phi=1$) with aqueous volume $\Omega^{0}$. The molar
concentration of OA vesicles in the extruded suspension is approximately

$C_{\text{ves}}^{\text{OA}}=\frac{RC_{0}}{N^{\text{OA}}}$ (17)

where $N^{\text{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]_{\text{eq}}^{\text{OA}}$, 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=\frac{\rho}{1-\rho}$) 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 ($\Phi=1$) with
aqueous volume $\Omega^{0}$, but now part of the bilayer consists of DOPA
phospholipid in molecular fraction $\rho$. The molar concentration of DOPA:OA
vesicles in the extruded suspension is approximately

$C_{\text{ves}}^{\text{DOPA:OA}}=\frac{C_{0}(1+g)}{N^{\text{DOPA:OA}}}$ (18)

where the aggregation number $N^{\text{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]_{\text{eq}}^{\text{DOPA:OA}}$, 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 $\Omega_{p}$ of each solution, creating a new mixed
system of volume $\Omega_{e}=2\Omega_{p}$, containing DOPA:OA vesicles in
number $N_{A}\Omega_{p}C_{\text{ves}}^{\text{DOPA:OA}}$ and OA vesicles in
number $N_{A}\Omega_{p}C_{\text{ves}}^{\text{OA}}$. The initial lipid monomer
concentration in the environment becomes
$\frac{1}{2}([L]_{\text{eq}}^{\text{DOPA:OA}}+[L]_{\text{eq}}^{\text{OA}})$.
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 $\Omega^{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 $\Phi=1$ in molar concentration

$C_{\text{ves}}^{\text{isotonic}}=\frac{C_{0}}{N^{\text{OA}}}$ (19)

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}\geq[B]+0.7$. This suspension is made of unilamellar
OA vesicles at $\Phi=1$ in molar concentration
$RC_{\text{ves}}^{\text{isotonic}}$, 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]_{\Delta}^{\text{max}}=0.16$M
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 $\Omega=\Omega^{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 $\Omega_{e}=2\Omega_{p}$ where isotonic vesicles number
$N_{A}\Omega_{p}C_{\text{ves}}^{\text{isotonic}}$, and the swelled vesicles
number an $R$ multiple of this. The lipid monomer concentration in this new,
larger environment is initially $[L]_{\text{eq}}^{\text{OA}}$.

#### 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 $\epsilon\approx 0.21$ for our model
pure oleate vesicles, where these vesicles burst through excessive osmotic
pressure when $\Phi<1-\epsilon$. Pure OA vesicles reached a minimum of
$\Phi=0.77$ in our phospholipid-driven competition simulations, and a minimum
of $\Phi=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 acid333General
note: The above procedures define a ‘concentration approach’ to mixing, where
two equal volumes are mixed, and the number of vesicles in the variable
population is controlled by increasing or decreasing vesicle concentration.
Another approach to mixing would be the ‘volume approach’ whereby the variable
population has a fixed vesicle concentration, but instead a variable volume
which controls the number of vesicles present. Volume mixing was found to
produce nearly equivalent outcomes, so only results following the
concentration mixing procedure are here reported..

##  III Results

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

Figure 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
($\rho=0.1$) vesicles are mixed 1:1 with pure OA vesicles (following our
mixing procedure with $\Omega_{p}=\Omega_{\text{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.
Figure 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<\rho_{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 $\rho_{0}^{\text{crit}}>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<\rho_{0}<0.25$, (ii) $0<\rho_{0}<0.5$, (iii) $0.25<\rho_{0}<0.75$ and
(iv) $0.3<\rho_{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\leq[B]_{\Delta}<0.16$ molar. Any vesicle starting at tension state
$\Phi_{0}^{\text{crit}}>0.802$ gains lipids from its competitors. See text for
full discussion.

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 ($\rho=0.1$) vesicles, when
this population is mixed $1:1$ with (i) pure OA vesicles (green lines), (ii)
DOPA:OA ($\rho=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 ($\rho=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 ($\rho=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
($\rho=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 phospholipid-driven competition in our
model is for vesicles stealing lipid to grow in surface and finish at high
$\Phi>1$ values (excess surface, flaccid), and for vesicles losing lipid to
suffer reduced surface, finishing at $\Phi<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 osmotically-driven competition can finish at a variety of surface
sizes. However, unlike phospholipid-driven competition, all OA vesicles will
finish with the same $\Phi<1$ value (equal osmotic stress). Here, surface
changes do not translate to final differences in $\Phi$, 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 present444The
initial lipid $L$ content of each individual vesicle is not explicitly part of
these boundary conditions. In fact, in our model, when total lipid $L_{t}$ is
fixed, initial vesicle surface sizes have no effect on the final equilibrium
of the system, only on the transient dynamics leading there.. 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 residual buffer up to 0.16M, vesicles with low
initial $\Phi$ values steal lipid from those with higher (less swelled) $\Phi$
values, with the tipping point between growing and shrinking at
$\Phi_{0}^{\text{crit}}=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
($\rho=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
$\rho^{\text{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.

Parameter | Description | Value | Unit
---|---|---|---
$S_{\mu}^{0}$ | Surface area of 100nm spherical vesicle | $3.142\times 10^{4}$ | [$nm^{2}$]
$\Omega^{0}$ | Volume of 100nm spherical vesicle | $5.236\times 10^{5}$ | [$nm^{3}$]
$[B]$ | Buffer concentration | $0.2$ | [$M$]
$k_{\text{out}}$ | OA lipid release constant | $7.6\times 10^{-2}$ | [$s^{-1}$]
$k_{\text{in}}$ | OA lipid uptake constant | $7.6\times 10^{3}$ | [$s{}^{-1}$$M^{-1}$$nm^{-2}$]
$\alpha_{L}$ | OA lipid head area(1) | 0.3 | [$nm^{2}$]
$\alpha_{P}$ | DOPA lipid head area(1) | 0.7 | [$nm^{2}$]
$\epsilon$ | OA vesicle burst tolerance | $0.21$ |
$[L]_{\text{eq}}^{\text{OA}}$ | 100nm OA vesicle, OA monomer equilibrium concentration | $6.667\times 10^{-5}$ | $[M]$
$[L]_{\text{eq}}^{\text{DOPA:OA}}$ | 100nm DOPA:OA vesicle, OA monomer equilibrium concentration | $5.294\times 10^{-5}$ | $[M]$
$N^{\text{OA}}$ | 100nm OA vesicle aggregation number | 209439 | total lipids
$N^{\text{DOPA:OA}}$ | 100nm DOPA:OA vesicle aggregation number | 184799 | total lipids
$\Omega_{\text{stoi}}$ | Competition volume unit for stoichiometric calculations | $3.478\times 10^{-13}$ | [$dm^{3}$]
$\Omega_{\text{dyn}}$ | Competition volume unit for dynamics simulations | $6.956\times 10^{-15}$ | [$dm^{3}$]
$C_{0}$ | Mix concentration unit | 0.001 | [$M$]
$[B]_{\Delta}^{\text{max}}$ | Residual buffer concentration inside maximally swelled OA vesicles(2) | $0.16$ | $[M]$
Table 1:  Vesicle competition model parameters. Figure 6:  Wider exploration
of three different vesicle competition scenarios. Relative surface growths of
two vesicle populations is explored in a broader context for three different
competition scenarios detailed by the key. A Phospholipid competition.
Population 1, a fixed population of vesicles with initial DOPA phospholipid
fraction $\rho_{0}^{pop1}=0.1$, is mixed $1:R$ with population 2, whose
vesicles have initial DOPA fraction $\rho_{0}^{pop2}$. B Osmotic competition.
Population 1, a fixed population of vesicles swelled by residual buffer
$[B]_{\Delta}=0.08M$, is mixed $1:R$ with population 2, whose vesicles are
swelled by residual buffer $[B]_{\Delta}^{pop2}$. C Phospholipid versus
osmotically-driven competition. Vesicles with initial DOPA fraction
$\rho_{0}^{pop1}$ are mixed $1:1$ with pure oleate vesicles swelled by
residual buffer $[B]_{\Delta}^{pop2}$. In all cases, for DOPA laden vesicles,
the direct effect is maximally enabled ($d=1$). Blue lines on plots highlight
when the relative surface growth is 1.

##  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)555 Likewise, the assumption we make with osmotically-driven
competition is that the residual buffer inside the vesicles permeates very
slowly through the bilayer membrane.. 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]_{\text{eq}}$, which is
trivially when all vesicles have the same number of membrane phospholipids
$P_{\mu}$. With no remaining asymmetries in $P_{\mu}$ 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_{\mu}$ 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 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
($\rho=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_{\mu}$ 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$, (10)) tends to _but
crucially never actually reaches_ the CVC concentration of pure oleic acid:

$\underset{L_{\mu}\rightarrow\infty}{lim}f=\frac{2k_{\text{out}}}{k_{\text{in}}\alpha_{L}}=[L]_{\text{eq}}^{OA}$
(20)

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]_{\text{eq}}^{OA}$. 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

$\textbf{u}({\Phi})=exp\left(\frac{1}{\Phi}-1\right)$ (21)

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

$\underset{L_{\mu}\rightarrow\infty}{lim}f=\frac{2k_{\text{out}}}{k_{\text{in}}\alpha_{L}}\cdot
exp(1)>[L]_{\text{eq}}^{OA}$ (22)

meaning that the DOPA:OA vesicles do not show the same continued growth as the
lipid monomer concentration in the environment rises toward
$[L]_{\text{eq}}^{OA}$. 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]_{\text{eq}}^{OA}$. 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é et al., (2007); Solé, (2009). 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.

## Acknowledgments

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 590-13), Spanish Ministry of Science (MINECO Grant FFI2011-25665),
COST Action CM 1304 (Emergence and Evolution of Complex Chemical Systems).
F.M. acknowledges support from MIUR (PRIN 2010/11 2010BJ23MN_003).

## References

  * Adamala and Szostak, (2013) Adamala, K. and Szostak, J. W. (2013).  Competition between model protocells driven by an encapsulated catalyst.  Nat Chem, 5:495–501.
  * Budin and Szostak, (2011) Budin, I. and Szostak, J. (2011).  Physical effects underlying the transition from primitive to modern cell membranes.  Proc Natl Acad Sci U S A, 108(13):5249–5254.
  * Cape et al., (2011) Cape, J. L., Monnard, P.-A., and Boncella, J. M. (2011).  Prebiotically relevant mixed fatty acid vesicles support anionic solute encapsulation and photochemically catalyzed trans-membrane charge transport.  Chem Sci, 2:661–671.
  * Chen et al., (2004) Chen, I. A., Roberts, R., and Szostak, J. W. (2004).  The emergence of competition between model protocells.  Science, 305:1474–1476.
  * Dyson, (1985) Dyson, F. (1985).  Origins of Life.  Cambridge University Press.
  * Gánti, (1975) Gánti, T. (1975).  Organization of chemical reactions into dividing and metabolising units: the chemotons.  Biosystems, 7:15–21.
  * Gillespie, (1976) Gillespie, D. T. (1976).  A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.  J Comput Phys, 22:403–434.
  * Lotka, (1925) Lotka, A. J. (1925).  Elements of Physical Biology.  Williams and Wilkins, Baltimore.
  * Luisi, (2006) Luisi, P. L. (2006).  The Emergence of Life: From Chemical Origins to Synthetic Biology.  Cambridge University Press.
  * Macía and Solé, (2007) Macía, J. and Solé, R. V. (2007).  Synthetic turing protocells: vesicle self-reproduction through symmetry-breaking instabilities.  Philos Trans R Soc Lond B Biol Sci, 362:1821–1829.
  * Mavelli and Ruiz-Mirazo, (2007) Mavelli, F. and Ruiz-Mirazo, K. (2007).  Stochastic simulations of minimal self-reproducing cellular systems.  Philos Trans R Soc Lond B Biol Sci, 362:1789–1802.
  * Mavelli and Ruiz-Mirazo, (2010) Mavelli, F. and Ruiz-Mirazo, K. (2010).  Environment: a computational platform to stochastically simulate reacting and self-reproducing lipid compartments.  Phys Biol, 7:036002.
  * Mavelli and Ruiz-Mirazo, (2013) Mavelli, F. and Ruiz-Mirazo, K. (2013).  Theoretical conditions for the stationary reproduction of model protocells.  Integr Biol, 5:324–341.
  * Mouritsen, (2005) Mouritsen, O. G. (2005).  Life - As a Matter of Fat: The Emerging Science of Lipidomics.  Springer-Verlag.
  * Munteanu et al., (2007) Munteanu, A., Attolini, C. S.-O., Rasmussen, S., Ziock, H., and Solé, R. V. (2007).  Generic darwinian selection in catalytic protocell assemblies.  Philos Trans R Soc Lond B Biol Sci, 362:1847–1855.
  * Rasmussen et al., (2009) Rasmussen, S., Bedau, M., Chen, L., Deamer, D., Krakauer, D., Packard, N., and Stadler, P., editors (2009).  Protocells: Bridging Non-Living and Living Matter.  MIT Press.
  * Ruiz-Mirazo et al., (2014) Ruiz-Mirazo, K., Briones, C., and de la Escosura, A. (2014).  Prebiotic systems chemistry: New perspectives for the origins of life.  Chem Rev, 114(1):285–366.
  * Ruiz-Mirazo et al., (2011) Ruiz-Mirazo, K., Piedrafita, G., Ciriaco, F., and Mavelli, F. (2011).  Stochastic Simulations of Mixed-Lipid Compartments: From Self-Assembling Vesicles to Self-Producing Protocells, volume 696 of Adv Exp Med Biol, pages 689–696.  Springer.
  * Segré et al., (2001) Segré, D., Ben-Eli, D., Deamer, D. W., and Lancet, D. (2001).  The lipid world.  Orig Life Evol Biosph, 31:119–145.
  * Smith and Szathmáry, (1995) Smith, J. M. and Szathmáry, E. (1995).  The Major Transitions in Evolution.  Oxford University Press.
  * Solé, (2009) Solé, R. V. (2009).  Evolution and self-assembly of protocells.  Int J Biochem Cell Biol, 41:274–284.
  * Solé et al., (2007) Solé, R. V., Munteanu, A., Rodriguez-Caso, C., and Macía, J. (2007).  Synthetic protocell biology: from reproduction to computation.  Philos Trans R Soc Lond B Biol Sci, 362:1727–1739.
  * Varela et al., (1974) Varela, F. G., Maturana, H. R., and Uribe, R. (1974).  Autopoiesis: the organization of living systems, its characterization and a model.  Biosystems, 5:187–196.