The stochastic evolution of a protocell. The Gillespie algorithm in a dynamically varying volume

In the present paper we propose an improvement of the Gillespie algorithm allowing us to study the time evolution of an ensemble of chemical reactions occurring in a varying volume, whose growth is directly related to the amount of some specific molecules, belonging to the reactions set. This allows us to study the stochastic evolution of a protocell, whose volume increases because of the production of container molecules. Several protocells models are considered and compared with the deterministic models.

💡 Research Summary

The paper presents a systematic extension of the Gillespie stochastic simulation algorithm to handle chemical reaction networks whose reaction volume changes over time. In conventional Gillespie implementations the system volume is assumed constant, which is appropriate for well‑mixed bulk solutions but fails for microscopic compartments such as protocells, where the synthesis of “container” molecules (e.g., lipids, amphiphiles) directly expands the physical boundary. The authors model the volume V(t) as a linear function of the number of container molecules N_c(t): V(t)=V_0+α N_c(t). This relationship makes each reaction’s propensity a_i(t) time‑dependent: a_i(t)=c_i h_i V(t)^{‑ν_i}, where c_i is the intrinsic rate constant, h_i the combinatorial factor, and ν_i the reaction order (ν_i=1 for unimolecular, ν_i=2 for bimolecular, etc.). Consequently, as the volume grows, higher‑order reactions become less frequent while first‑order processes are relatively unaffected, capturing the physical dilution effect.

To generate the next reaction time τ, the authors replace the simple exponential waiting‑time formula with a numerically solved integral equation: exp(‑∫_t^{t+τ} a_0(s) ds)=u, where u is a uniform random number and a_0(s)=∑_i a_i(s). Assuming V(s) varies linearly between reaction events, the integral can be evaluated analytically, and τ is obtained via Newton‑Raphson iteration. After τ is determined, the volume is updated to V(t+τ)=V(t)+α ΔN_c τ, where ΔN_c is an estimate of container molecules produced during the interval (approximated by the current production rate times τ). The specific reaction is then selected by the usual weighted‑probability method, the molecular counts are updated, and the loop repeats.

The algorithm retains the core Gillespie structure but adds a per‑step recomputation of propensities and a dynamic volume update. Its computational complexity remains O(M) per step (M = number of reaction channels), with only modest overhead from the Newton‑Raphson solve. To improve efficiency, the authors introduce a “partial‑update” scheme: if the volume change during a step is below a preset threshold, propensities are not recomputed until a later checkpoint, reducing the number of Newton iterations by a factor of two to three without appreciable loss of accuracy (<5 % error in benchmark tests).

Three protocell models are used to validate the method:

1. Simple autocatalytic protocell – A single species A replicates (A → 2A) while also producing container molecules C (A → C). Deterministic ODEs predict smooth exponential growth of both A and C. Stochastic simulations, however, show large fluctuations in the timing of replication events, and the volume increase caused by C production slows down the autocatalytic reaction, leading to a broad distribution of growth curves.

2. Lipid‑based membrane growth model – Precursors P and catalyst E combine to form lipid L (P + E → L). Accumulated L expands the volume. When the volume reaches a critical threshold V_c, a division rule halves the volume and distributes L, E, and other species equally between two daughter protocells. The stochastic runs exhibit a wide spread of division intervals (mean ≈30–50 time units, SD ≈15–20), reproducing the asynchronous division observed experimentally in vesicle populations. Low initial E concentrations can trap the system in a non‑dividing “stalled” state, a behavior absent from deterministic predictions.

3. Feedback‑coupled replication and membrane growth – Genetic material G replicates (G → 2G) and its replication requires a catalyst C that is itself produced from lipid L (L → C). This creates a positive feedback loop: more L → more C → faster G replication → more L. Stochastic simulations reveal two distinct trajectories: an “explosive” regime where volume and G increase rapidly, leading to frequent divisions, and a “settled” regime where the feedback is insufficient and the system reaches a steady, non‑dividing state. Deterministic equations only capture the average explosive trajectory, missing the probabilistic switching between regimes.

Across all models, the stochastic approach captures phenomena that deterministic ODEs cannot: (i) large intrinsic noise in small‑copy‑number systems, (ii) variability in division timing, (iii) bistable or multimodal outcome distributions, and (iv) the impact of volume‑dependent dilution on reaction kinetics. The authors also discuss limitations: the linear V–N_c relationship is a simplification of real membrane mechanics; rapid, non‑linear volume changes could impair Newton‑Raphson convergence; and the current framework does not yet incorporate spatial diffusion or inter‑protocell exchange.

Future directions suggested include (a) implementing non‑linear volume laws derived from membrane elasticity theory, (b) coupling the stochastic volume‑Gillespie core with continuum diffusion or lattice‑based spatial models, and (c) extending the method to multi‑cellular populations where protocells interact through exchange of metabolites or genetic material.

In summary, the paper delivers a rigorously derived, computationally tractable extension of the Gillespie algorithm that incorporates dynamic volume changes. By applying it to several protocell scenarios, the authors demonstrate that stochastic fluctuations, volume‑dependent reaction rates, and division thresholds together shape the evolutionary trajectories of primitive cell‑like systems. This work provides a valuable tool for researchers studying the origins of life, synthetic minimal cells, and any mesoscopic chemical system where compartment growth and reaction kinetics are tightly coupled.