# Continuous Game of Life: a minimal model for the emergence of cell division and motility

Alexandre Guillet, Frank Jülicher

<figure><img src="image_1.png" /></figure>

## To cite this version:

Alexandre Guillet, Frank Jülicher. Continuous Game of Life: a minimal model for the emergence of cell division and motility. Proceedings of the Artificial Life Conference 2025, Oct 2025, Kyoto, Japan. hal-05324712

HAL Id: hal-05324712

https://hal.science/hal-05324712v1

Submitted on 21 Oct 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Continuous Game of Life:
a minimal model for the emergence of cell division and motility

Alexandre Guillet and Frank Jülicher

Max Planck Institute for the Physics of Complex Systems, Dresden, Germany
aguillet@pks.mpg.de

## Abstract

Conway's Game of Life shows that simple rules can generate a rich diversity of emerging structures. More recently, this cellular automaton has been translated to continuous space by (Rafler, 2011) in a simulation called SmoothLife. Related to its isotropic rule, the beauty of patterns generated by this continuous Game of Life has attracted the attention of a growing community at the intersection of science and computer art. We study a minimal variant of the continuous Game of Life that generates cell-like patterns capable of self-replicating, gliding and disappearing. Variations of the emerging dynamic morphology of these unit patterns are described from the perspective of bifurcations, instabilities and dynamical phase transition. Its mapping onto a large reaction-diffusion system is interpreted in terms of homeostatic concentrations of morphogens, maintained by a cascade of fast auxiliary reactions. Introducing a global conservation law, which constrains resource availability and consumption, is shown to induce self-organized criticality. Resemblance to biological processes such as division, motility, and death—together with a concise formulation and numerical implementation—makes it an appealing toy model for investigating the emergence of life-like patterns.

Code available at: codeberg.org/A-Guillet/cGoL

## Introduction

Since its inception in 1970, the cellular automaton *Game of Life* (GoL) defined by Conway (Gardner, 1970), famous for its complexity out of simple rules on the 2-dimensional grid, has fostered a strong interest and numerous developments. After the pioneering efforts in porting the Game of Life to the continuous space (MacLennan, 1990), in particular its glider pattern, via a scaling limit of large (inner and outer) discrete neighbourhoods (Evans, 2001, 2003; Pivato, 2007), the quest for more realistic Life-like automata has slowly faded away from academia.

In contrast, the exploration of the Game of Life became active in online communities of enthusiasts<fnref n="1" /> supported by the popularization of computers. This activity has led to the construction of increasingly complex structures within the rules of the original GoL model, such as a functioning

computer and an abstract self-replicating machine, illustrating the consequences of its Turing-completeness (Rendell, 2002).

The continued interest for the GoL and its variants has resulted in the first successful implementation of a space-continuous isotropic rule, called *SmoothLife* (Rafler, 2011), based on disk and annulus-shaped integration domains in place of neighbour counting regions. Simulations of gliders with better continuity in time have quickly followed, based on the explicit Euler method modified with a clamp regularization (Rafler, 2012; Hutton et al., 2015). This breakthrough has set the basis for a renewed exploration of Life-like continuous systems inspired from cellular automata but cleared from discrete states and the 4-fold symmetry of the grid. We call them *continuous Games of Life* (cGoL).

Among *SmoothLife* variations, several conceptual simplifications occurred, while maintaining a rich phenomenology. Most relevant for this article, the little known *SmootherLife* refinement from CornusAmmonis unifies the inner and outer integrals as convolutions with respect to a small kernel and a large kernel with the same Gaussian shape (CornusAmmonis, 2017; Hutton et al., 2015). This model generates clear self-replicating unit patterns made of a nucleus and a shell that resemble dividing (and sometimes dying) biological cells. This striking observation has motivated our investigation.

Another variation consists in reducing the nonlinearity of the GoL's survival rule, from its original bivariate form to a univariate function of the convolution integral. Called *Lenia*, this simplification comes at the price of a more sophisticated shape for the unique isotropic convolution kernel (Chan, 2019). The resulting simulation can generate single gliders (or oscillators) with surprisingly diverse morphologies, leading to the proposal of a naturalistic classification for artificial lifeforms. This exploration is pursued in (Chan, 2020), where the extended model returns to a multivariate survival function (restricted to a sum of univariate ones), along with other higher-dimensional extensions.

The *Lenia* model has generated a strong interest, uncovering certain caveats and fostering developments. In spite

<footnote n="1">¹See conwaylife.com and the Golly software.</footnote>
of the small time steps, discreteness in time still plays an important role in the stability of the patterns: *Lenia's* gliders tend to be destabilized and vanish in the continuous-time limit (Davis and Bongard, 2022; Kojima and Ikegami, 2023). When reformulating the *Lenia* model without the clamp regularization step, gliders can still be found in the continuous-time limit although with differing morphologies (Kawaguchi et al., 2021). The mathematical implications of these two modelling approaches are discussed in (Calcaterra and Boldt, 2023). The resulting integro-differential equation can be interpreted in terms of a particular reaction-diffusion system (Kojima and Ikegami, 2023): the slow field reacts to a combination of many fast diffusing and decaying species that it generates.

Furthermore, *Lenia's* gliders are metastable local structures, coexisting with stable global solutions: they easily vanish or explode into a space-filling pattern upon interaction. In order to address this issue, an environmental retroaction from resource rarefaction has been introduced in (Suzuki et al., 2023). The *flow-Lenia* model (Plantec et al., 2023) proposes a complete reformulation of the extended *Lenia* model (Chan, 2020) in terms of fluxes, which enforces a strict conservation law.

In this article, we integrate all these approaches into a minimal model that offers a continuous formulation of the discrete Game of Life while producing clear self-replicating and motile elementary patterns.

We first introduce the cGoL model: after a generic formulation of the rule of the game, we specify a parametrization with a minimal number of parameters. The reference choice of 6 + 1 fine-tuned values gives rise to cell-like patterns in two dimensions that either divide, or glide and oscillate, while nearby parameter values yield a broader diversity of behaviours and morphologies. These variations are discussed in terms of linear stability, bifurcations, shape symmetry and instabilities, and dynamical phase transition. We then examine physical aspects of the model. Its interpretation as a reaction-diffusion system can account for the Gaussian shape of the convolution kernels through a cascade of auxiliary reactions. Finally, we introduce a global conservation law that constrains resource availability and consumption via feedback mechanisms, leading to self-organized criticality. The article concludes with perspectives and implications of this formulation.

# Definition of a continuous Game of Life

## Generic model

Let $L(x, t)$ be the field of a continuous Game of Life (cGoL), in Euclidean space $x \in \mathbb{R}^d$ and time $t > 0$, taking scalar values in $\mathbb{R}$. The continuous dynamics is determined by a

nonlinear partial integro-differential equation of the form

$$ \partial_t L = \Gamma(\Phi_k * L) - L \quad (1) $$

$$ \Phi_k * L(x, t) = \int \Phi_k(x - y)L(y, t) dy \quad . \quad (2) $$

Here, $\Gamma$ is a nonlinear function of multiple variables indexed by $k = 1, 2, \dots$, extending the univariate (*Lenia*) case formulated in Kawaguchi et al. (2021). The variables $\Phi_k * L$ are auxiliary fields that represent neighbourhood population states, related to the primary field $L$ by a linear (non-local) operation in space. The integral Eq. (2) serves as the continuous analogue to counting discrete neighbours. It represents a local average of $L$ around any spatial coordinate $x$, which is reflected in its convolution structure.

Convolution kernels $\Phi_k$ correspond to different types of neighbourhood, typically a small one ($k = 1$) and a large one ($k = 2$) about 3 times wider. Isotropic kernels, $\Phi_k(x) = \Phi_k(|x|)$, are a natural choice to restore the symmetry lacking in the discrete model. A constant kernel amounts to a simple integration, and will be introduced through a global conservation law ($k = 3$).

The resulting growth term $\Gamma(\Phi_k * L)$ is called the target field, moving as $L$ relaxes towards it. The timescale of this relaxation is taken as the time unit in Eq. (1). Defining the target of the continuous Game of Life amounts to specifying the discrete rule in the original Game of Life.

## Specific model

In this article, we remain close to the original rule of the *Game of Life* in two ways: by starting with a bivariate target function (survival rule), as in *SmoothLife*, and by using the original scaling and the same shape for convolution kernels (neighbourhoods) as in *SmootherLife*.

Instead of discrete squares (of widths 1 and 3), the Gaussian shape describes neighbourhoods as isotropic and smooth kernels:

$$ \Phi_1(x) = e^{-\pi|x|^2} \quad , \quad \Phi_2(x) = \lambda^{-d}\Phi_1(x/\lambda) \quad . \quad (3) $$

Therefore, the resulting auxiliary fields, which we denote $M = \Phi_1 * L$ and $N = \Phi_2 * L$, are simple Gaussian blurs of the primary field, as illustrated in Fig. 1 (top). We take the smallest scale as the unit length and set the larger one to the original scaling factor $\lambda = 3$.

By evaluating the target function $\Gamma$ for pairs of values $(M, N)$ at every spatial coordinate, one obtains the target field $\Gamma(M, N)$ that drives the dynamics Eq. (1) of the primary field $L$, as illustrated in Fig. 1 (bottom). The target is chosen to take values in $[0, 1]$, so that $L$ is eventually constrained in this range (after a short transient when initialized or perturbed outside it). Originally defined as a Boolean table, $\Gamma$ describes a survival range in the population of the largest neighbourhood $N$, which varies with the state of the
<figure><img src="image_3.png" /><figcaption>Figure 1: Primary and auxiliary fields: L, M = $\Phi_1 * L$ and $N = \Phi_2 * L$ (top, from left to right). Target function $\Gamma$ (bottom left, same gray scale as for L, M, N). Special uniform states are marked in red: post-bifurcation unstable and stable (empty and full triangles) steady states, and pre-bifurcation non-steady optimum (empty circle). Spatial distribution of pairs (M, N), where $\Gamma$ gets evaluated (bottom centre, pixel count). Growth or decay rate of the field L, Eq. (1) (bottom right).</figcaption></figure>

smallest neighbourhood M. Building upon previous definitions that parametrize the kink of the discrete rule, we retain the form represented in Figs. 1 and 2,

$$ \Gamma(M, N) = S'\left(\frac{N-N_c(M)}{\delta N_c(M)}\right) \quad (4) $$

$$ N_c(M) = N_0 + (N_1 - N_0)S\left(\frac{M-M_c}{\delta M_c}\right) $$

$$ \delta N_c(M) = \delta N_0 + (\delta N_1 - \delta N_0)S\left(\frac{M-M_c}{\delta M_c}\right) $$

$$ S(m) = (1 + e^{-4m})^{-1}, \quad (5) $$

where the derivative $S'$ of the sigmoid is a bell shape: $S'(n) = 4S(n)(1 - S(n)) = \text{sech}(2n)^2$. The target parameters are specific values of the auxiliary fields M and N denoted with an index, see Eq. (7). This definition of the survival range effectively diminishes parametric degrees of freedom by one as compared to (Rafler, 2011; CornusAmmonis, 2017). The locality of the growth region appears crucial, so that other fast decaying bell shapes such as the Gaussian function yield comparable results.

The equivalent form

$$ \Gamma(M, N) = S'(n_0 + (n_1 - n_0)S(\tilde{m})) \quad (6) $$

$$ (\tilde{m}, n_0, n_1) = \left( \frac{M - \tilde{M}_c}{\delta M_c}, \frac{N - N_0}{\delta N_0}, \frac{N - N_1}{\delta N_1} \right) $$

arises from properties of the sigmoid function S, Eq. (5). Note that the first parameter differs as compared to Eq. (4): $\tilde{M}_c = M_c - \frac{\delta M_c}{4} \log \frac{\delta N_1}{\delta N_0}$.

The rule of the Game is fully specified by the shapes $\Gamma$ (based on the sigmoid S) and $\Phi$ (Gaussian), and their parameters ($M_c$, $\delta M_c$, $N_0$, $\delta N_0$, $N_1$, $\delta N_1$) and $\lambda$. Given these parametrizations of the target function and kernels, the parameter space thus has 6 + 1 dimensions (to which one can add the dimensionality of space d). In this article, we use Eqs. (4–5) in space with d = 2 dimensions, original scaling $\lambda = 3$, and reference target parameters p:

$$ \begin{aligned} \mathbf{p} &= (M_c, \delta M_c, N_0, \delta N_0, N_1, \delta N_1) \\ &= (0.50, 0.10, 0.23, 0.015, 0.35, 0.26) . \end{aligned} \quad (7) $$
These parameters have been fine-tuned (and rounded to the closest values with two significant digits) so that the elementary cell-like pattern is both able to divide several times and to glide while oscillating, as illustrated in Fig. 2 and described next.

## Emergence of dynamic patterns

### Uniform states and bifurcations

Some aspects of this system can be discussed from the scalar case $d = 0$, which corresponds to uniform solutions for $d > 0$. In this case, the convolution kernel is just a coefficient, equal to one since normalized in Eq. (2): $L = M = N$. Near the reference parameters Eq. (7), the system has 3 uniform steady states (fixed points), marked with triangles in Fig. 1. These solutions of the equation $\Gamma(L_*, L_*) = L_*$ are practically undistinguishable from the solutions of $S'(\frac{L_* - N_0}{\delta N_0}) = L_*$ as they occur in the lower region of the target before the transition ($N_0 < M_c - 2\delta M_c$).

Since the survival range in this region is very thin, $\delta N_0 \ll N_0$, the lowest (empty) uniform steady state is practically vanishing, and it is stable (Game over). The two other ones lie in the vicinity of $N_0 \pm \delta N_0$, hence they are very close to one another and related to the shell of the cell pattern. This pair of unstable and stable fixed points may be understood as tuning the proximity with a saddle-node bifurcation. An empty state with exact zero value, or a crossing of the bifurcation would both require a reparametrization of $\Gamma$.

For different target parameters, another saddle-node bifurcation can occur. The local maximum of $\Gamma(L, L) - L$ corresponds to the non-steady uniform state closest to the bifurcation. This special value, marked with a circle in Fig. 1, depends on the four other parameters, $M_c, \delta M_c, N_1, \delta N_1$, and relates to the nucleus of the cell pattern.

The parrot-shaped distribution of values $(M, N)$ shown in Fig. 1 (bottom centre), which corresponds to a non-uniform non-steady state $L$, has visible traces of these special uniform values as cusps or recess near their location. The rich morphogenesis generated by this nonlinear system exhibits numerous other bifurcations, whose analytical exploration proves challenging. Consequently, we opt for an empirical approach to investigate these phenomena.

### Emerging morphologies via instabilities

Near the reference parameters Eq. (7), the continuous Game of Life has a diversity of steady and non-steady local non-uniform solutions, which are all variants of a recognizable elementary pattern. Activated regions $L$ form a nucleus and a shell, which correspond respectively to the thick and the thin regions of the survival range in the target $\Gamma$ (for high and low $M$). This correspondence is colour-coded in Fig. 2. The inner and outer empty areas correspond to the regions at the right and at the left of the survival range respectively. In analogy with biology, we call this elementary pattern a *cell*, due to its capabilities of self-replication, motility and death.

<figure><img src="image_4.png" /><figcaption>Figure 2: Cell dynamics at and near reference parameters Eq. (7), with colour coding of corresponding regions of the target function $\Gamma$ (A) and the field $L$ (B-E). Reference behaviours: dividing cell (B) or oscillating glider (C), depending on its symmetry axes (top). Time is incremented from top to bottom of 10 (B) and 6 (C) time units. (D-E) Variants: (steady) isotropic egg (D) and (clockwise) turning glider (E) for $N_1$ respectively increased and decreased by 0.03.</figcaption></figure>

In unfavourable regions of the parameter space where the unit pattern can still be initialized<fnref n="2" />, the cell remains in the form of a dormant *egg*, which is an isotropic steady solution shown in Fig. 2 (D), or may die (vanish). Bistability can already be observed, with two possible radii of the egg at certain parameter values.

For more favourable parameters, the egg grows and can undergo various shape instabilities, triggered by environmental perturbations (such as interaction with neighbour cell patterns). The first symmetry breaking leads to the polarization of the egg with a long and a short symmetry axes. The elongation of the cell can lead to its *division*, while keeping its two symmetry axes, see Fig. 2 (B). Some parameters lead to *failed division*, in which case the cells die shortly after the division of the nucleus. If isotropy is not broken, the egg may grow into a *bubble*.

The short symmetry axis can turn unstable, leading to a second symmetry breaking. The result is a *glider*: a bilateral

<footnote n="2">²We use the seed $L(0, \mathbf{x}) = e^{-\pi \mathbf{x}^T \Sigma^{-1} \mathbf{x}}$ with a slight elongation to facilitate the first division by breaking isotropy.</footnote>
cell morphology associated to motion. The glider can have a variety of behaviours and shape instabilities. In the reference case illustrated in Fig. 2 (B, C), the short axis is metastable so that dividing cells can be destabilized into a glider. Conversely, interacting gliders can stop and divide. This glider is notable for generating shape oscillations during motion. When amplified, this oscillation leads back to cell division, with intermediary cases for which division fails asymmetrically. For other parameters, oscillations can be damped and result in a glider with a steady shape.

The long symmetry axis can also turn unstable, leading to a last symmetry breaking. It typically yields a clockwise or anti-clockwise *turning glider*, such as Fig. 2 (E), with various radii for its circular trajectory. A *tripod*, steady cell with order 3 symmetry, can also be found, and observed to switch its morphology to different gliders upon perturbation.

The shape instability of an elongating cell may also be absent or insufficient to produce a successful division. This situation yields various extended patterns. If the division of the nucleus is not followed by the scission of the shell, an intermediary body can be formed in the centre, which may turn into a full nucleus. Otherwise, the formation of segmented *chains* or non-segmented *filaments* is initiated and can percolate. A nucleus that elongates but has difficulty to divide leads to a *snake*, with growing or shrinking ends. An ever-growing snake can fill space and generate various mazes and Turing patterns.

The phenomenology close to the reference cGoL system is unexpectedly diverse, in spite of the simple morphology of the elementary cell as compared to *Lenia's* gliders.

## Dynamical phase transition

The rich diversity of shape instabilities and examples of bistable morphologies in the direct vicinity of the reference parameters is reminiscent of criticality. We introduce a prototypical order parameter to simplify our understanding of the 6 parametric degrees of freedom in the target function, and we attempt a higher-dimensional visualization.

Consider the scaling of the target function, $r\Gamma$, with a constant coefficient $r$ (close to 1): $r > 1$ enhances the growth rate whereas $r < 1$ reduces it, according to Eq. (1). As a result, the field $L$ takes values in the interval $[0, r]$. Since the scaling coefficient $r$ is constant, this modulation of growth amounts to scaling all auxiliary fields instead, by redefinition of $L$ inside $[0, 1]$. In turn, this is equivalent to rescaling all target parameters, given the parametrization Eq. (4) of $\Gamma$. One obtains a simple interpretation of a perturbation along a particular direction of the parameter space:

$$r\Gamma(M, N; \mathbf{p}) \sim \Gamma(rM, rN; \mathbf{p}) = \Gamma(M, N; \mathbf{p}/r) \quad (8)$$

where $\sim$ denotes the equivalence of the resulting dynamics up to a rescaling of $L$. As expected, the scaling $\mathbf{p}/r$ Eq. (7) promotes or inhibits growth, for $r > 1$ and $r < 1$

respectively, interpolating across many of the previously described morphologies between two coarse phases of the system: dense or dilute.

<figure><img src="image_5.png" /><figcaption>Figure 3: Dynamic states of the cGoL model in the presence of a spatial variation of the target parameters $\mathbf{p} = \mathbf{p}(\mathbf{x})$. For each panel, a subset of the parameters (specified on the left) is increased linearly from left to right around its reference value Eq. (7) (centre). The top panel corresponds to Eq. (8) with a variable scaling $r = r(\mathbf{x})$ (top axis).</figcaption></figure>

We further attempt to visualize the role of individual target parameters or subsets of them on the morphogenesis by simulating the cGoL with a spatial parameter scans $\mathbf{p} = \mathbf{p}(\mathbf{x})$ along the horizontal direction. The resulting overview Fig. 3 maps adjacent patterns, such as mazes, filaments and uniform states, around the reference cell and reveals the sensitivity of the tuning with respect to the different target parameters.

Represented in the top panel of Fig. 3, the scaling Eq. (8) of all target parameters is dominated by the effect of scaling the auxiliary field $M$, in particular the first target parameter $M_c = \frac{1}{2}$. The isotropic egg and turning glider variants from Fig. 2 (D, E), with values of $N_1$ differing of nearly $\pm 9\%$, can be located on the map Fig. 3 in the second panel from the bottom. We can notice the propensity of this turning glider to collectively assemble in chains.

This visualization technique suggests that the reference cell morphology has been selected close to a dynamical phase transition, as will be quantified in a later section.
# Interpretation as a reaction-diffusion system
## Homeostatic concentrations of morphogens

As exposed in (Kojima and Ikegami, 2023) in the case of Le-
nia, a cGoL model can be mapped to an asymptotic regime
of certain reaction-diffusion systems

$$
\tau_k \partial_t M_k = \sigma_k^2 \nabla^2 M_k + f_k(M_{k'}) \quad (9)
$$

in which a primary field $M_0 = L$ does not diffuse ($\sigma_0 = 0$)
and reacts to some morphogens, among multiple auxiliary
species of concentration fields $M_k$ for $k \ge 1$, that diffuse
and react much faster: $\tau_k \ll \tau_0$. Extended models such
as (Chan, 2020; Plantec et al., 2023) consider multiple pri-
mary (slow) fields.

The only nonlinear reaction is $f_0 = \Gamma - L$, Eq. (1), in the growth rate $\Gamma$. Previous formulations based on *SmoothLife* also include $L$ in the nonlinearity, by modelling its saturation with a clamp function in the interval $[0, 1]$. The survival range encoded in $\Gamma$ represents the target morphogen concentrations, maintained by the dynamics, and mapped in the cell in Fig. 2. Therefore, target parameters Eq. (7) in the nonlinearity directly specify the homeostatic concentrations of morphogens in the cells.

The convolution form Eq. (2) for morphogen fields results
from the linear auxiliary reactions $f_k(M_{k'})$ for $k \ge 1$ and
multiple $k' \ge 0$, in which diffusion, decay and conversion
terms are allowed. Each fast-relaxed auxiliary field can be
solved in terms of the slow field $L$ and a certain Green func-
tion (its kernel): $M_k = \Phi_k * L$. The linearity may be under-
stood as describing small perturbations near homeostasis.

Cascade of fast auxiliary reactions

The Gaussian shape of the kernels Eq. (3) is the result of
a large cascade of fast auxiliary reactions $k = 1, 2, \dots$, as
captured by the simple linear reaction network:

$$
\tau \partial_t M_k = \sigma^2 \nabla^2 M_k - M_k + M_{k-1} . \quad (10)
$$

It describes a cascade of diffusing and decaying species tak-
ing its source in the cell’s shell and nucleus structures in
$M_0 = L$. Using the separation of timescales $\tau \ll 1$, relaxed
states are formally expressed with a linear space operator
$M_k = [1 - \sigma^2 \nabla^2]^{-k} L = \phi_k * L$. The corresponding ker-
nel shapes $\phi_k = \phi_1 * \cdots * \phi_1$ are obtained from Fourier
methods, transforming the Laplace operator $\nabla^2$ in terms of
the spatial frequency (wave vector), $-(2\pi|\mathbf{q}|)^2$. The inverse
Fourier transform yields their functional form

$$
\begin{align}
\phi_k(\mathbf{x}) &= \int \frac{e^{i2\pi\mathbf{x}\cdot\mathbf{q}}}{(1 + (2\pi\sigma|\mathbf{q}|)^2)^k} \mathrm{d}\mathbf{q} \nonumber \\
&= \frac{(|\mathbf{x}|/\sigma)^k K_{\frac{d}{2}-k}(|\mathbf{x}|/\sigma)}{2^{k-1}(2\pi\sigma|\mathbf{x}|)^{\frac{d}{2}}} \tag{11}
\end{align}
$$

involving the modified Bessel function of the second kind,
with exponential tail $K_\nu(x) \sim \sqrt{\frac{\pi}{2x}} e^{-x}$. This kernel cor-
responds to the steady spatial profile of the $k^{\text{th}}$ species in

the cascade, decaying and diffusing from a point source.
This Green function also assumes the absence of boundaries:
$x \in \mathbb{R}^d$. In the case of a short cascade, $k \le \frac{d+1}{2}$, a singular-
ity (or a cusp on the bound) is present at the origin.

The Gaussian shape arises for large cascades, $k \gg 1$, as
the consequence of the central limit theorem, summarized at
the level of the space operator as $[1 - \sigma^2\nabla^2]^{-k} \sim e^{k\sigma^2\nabla^2}$.
The Gaussian kernel with unit scale in Eq. (3) is matched
in the scaling limit $\sigma^2 \sim (4\pi k)^{-1}$. Large and small kernels
with scale ratio $\lambda = 3$ are idealized as morphogens at indices
$\lambda^2 k = 9k$ and $k \gg 1$ in the cascade Eq. (10).

Without cascade, the simplest Bessel shape $\phi_1$ represent
morphogens with different diffusion coefficients, directly
sourced from $L$ without intermediary reaction. In this case
gliding does not occur and divisions are incomplete, with
the systematic formation of an intermediary body. A small
cascade with $M = \phi_1 * L$ and $N = \phi_9 * L$ (or as small as
$\phi_4$) is sufficient for complete divisions.

In line with and extending the conclusions of (Kojima and Ikegami, 2023), the cGoL model can be viewed as a coarse-grained description of a large reaction-diffusion system, emphasizing the nonlinear response of the slow field to a few diffusing morphogens. Its distinctive feature is the non-standard nonlinearity $\Gamma$, inherited from the GoL survival rule, which directly sets the homeostatic morphogen concentration ranges. The Gaussian shape of the kernels follows from a linear cascade model of fast auxiliary reactions.

Limiting resource to self-organize at criticality
Introducing a global conservation law

In the reference cGoL system, cell division leads to un-
bounded growth, constrained only by the simulated volume.
To address this unphysical behaviour, we enforce a global
mass conservation on the system and ultimately show that it
gives rise to self-organized criticality.

Let us assume that the growth of the cells, in the slow
primary field $L$, consumes a finite resource, which only re-
plenishes when $L$ decays. A global conservation law can
be enforced between available resource $\bar{R}$ and consumed re-
source $\int L \, dx$ that constitutes the mass of the system<fnref n="3" />:

$$
R = \bar{R} + \int L \, dx \quad , \quad \partial_t R = 0 \quad . \quad (12)
$$

Drawing inspiration from the resource abundance field
in (Suzuki et al., 2023), a limitation of resource can be in-
troduced in the cGoL model via an availability coefficient

$$
r(t) = \frac{\bar{R}(t)}{R} \qquad (13)
$$

at most 1 for abundant resource (as previously), and vanish-
ing when the resource is fully consumed. A negative retroac-
tion on $L$ from decreasing globally available resource $\bar{R}(t)$

<footnote n="3">³The vanishingly small value of the empty state has been ne-
glected here, and should be subtracted from $L$ otherwise.</footnote>
can be achieved by inserting $r(t) \le 1$ in the target, as in Eq. (8).

The temporal modulation $r(t)$ breaks the equivalence in Eq. (8), yielding two distinct implementations: feedback type I (left-hand side) limits growth by decreasing the target $r(t)\Gamma$, whereas feedback type II (right-hand side) decreases auxiliary fields, $r(t)M$ and $r(t)N$, or equivalently increases target parameters, $p/r(t)$. Feedback of type II can be visualized as pushing the reference system towards the right in the top panel of Fig. 3, as the available resource is consumed. Both types I and II introduce the mass $\int L \, dx = 1*L$ as a third variable (with constant kernel) in the nonlinearity $\Gamma$.

By comparing the perturbation of different parameters in Fig. 3, one realizes that the resource rarefaction effect is dominated by the increase of the first parameter $M_c$. Therefore, we further devise feedback type III that preserves the bivariate formulation of the target function $\Gamma(M, N)$, but instead perturbs the kernel of the first auxiliary field

$$M = (\Phi_1 - \frac{M_c}{R}) * L, \quad (14)$$

where the constant negative term in the kernel yields a simple integration of $L$.

We then quantify the effect of coupling the field $L$ to the globally available resource $\bar{R}(t)$ from the conservation law Eq. (12), for each feedback type, respectively introduced in the dynamics Eq. (1) via $r(t)$ as in the left (I) and right-hand (II) sides of Eq. (8), or via Eq. (14) (III). Among many other possibilities, the retroaction mechanisms II and III can be directly interpreted as retuning, dynamically and autonomously, all reference parameters Eq. (7) or just the first one, with coupling strength $R^{-1}$.

## Phase transition from volume to resource limitation

Given a finite resource $R$, but also a necessarily finite spatial domain $\Omega \subset \mathbb{R}^d$ of volume $|\Omega|$, we count cells (at equilibrium$^4$) by thresholding the field $M$ at the transition value $M_c = \frac{1}{2}$ and by counting connected subsets of $\Omega$ satisfying $M > \bar{M}_c$ (distinct nuclei). This number may still fluctuate at long time, so that we estimate its mean value $\nu$ from both a temporal and an ensemble average across many simulations that have reached equilibrium. Similarly, the field $L$ relaxes to an equilibrium mean value $\mu$, estimated by averaging the density $\int L \, dx/|\Omega|$ in the same way. The mass per cell $m = \frac{\mu|\Omega|}{\nu}$ is then deduced.

Their dependence in the resource $R$ and volume $|\Omega|$ is summarized in Fig. 4: the number of cells $\nu$ (top left) is extensive with the limited resource $R$ provided a sufficient volume, or extensive with the limited volume $|\Omega|$ provided a sufficient resource. The discreteness of the elementary pattern is visible as steps for low cell count $\nu = 1, 2, 4, \dots$. The cross-over between the resource-limited and volume-limited regimes occurs around a specific value $\rho_*$ of the resource

<figure><img src="image_6.png" /><figcaption>Figure 4: Equilibrium states as a function of resource $R$ and volume $|\Omega|$, in terms of the number of cells $\nu$ (top left), the mean field $\mu$ (top right), and the mass per cell $m$ (bottom left). Lines encode different volumes (colours) and resource feedback types: target scaling (I), morphogen scaling (II), and kernel shift (III)</figcaption></figure>

density $R/|\Omega|$, as shown for the mean field $\mu$ in Fig. 4 (top right), well approximated as

$$\mu \approx \frac{\mu_{\infty}}{\sqrt{(\rho_{*} \frac{|\Omega|}{R})^2 + 1}} \quad (15)$$

where the asymptote for abundant resource $\mu_\infty = 0.131 \pm 0.001$ is a property of the reference cGoL model ($r=1$), whereas the estimates $\rho_*^{(\text{I, II})} = 10.5 \pm 0.5$ and $\rho_*^{(\text{III})} = 13.5 \pm 0.5$ are specific to the feedback type.

Irrespective of the approximation Eq. (15), the equality $\mu^{(\text{I})} = \mu^{(\text{II})} \Leftrightarrow r^{(\text{I})} = r^{(\text{II})}$ is observed to be robust in all regimes. Since feedback types I and II share the same long term resource abundance coefficient $r = r^{(\text{I, II})} = 1 - \mu^{(\text{I, II})} \frac{|\Omega|}{R}$, the same cell morphology is expected from Eq. (8), up to a rescaling of $L$. This rescaling explains the mismatch of the mass per cell at limited resource ($r < 1$): $m^{(\text{II})} = m^{(\text{I})}/r$, visible in Fig. 4 (bottom left) and rescaled in the inset. As a consequence, the number of cells also differs: $\nu^{(\text{II})} = r\nu^{(\text{I})}$.

The mass per cell varies from about $m_*^{(\text{II})} = m_*^{(\text{I})}/r_* = 6.10 \pm 0.05$ and $m_*^{(\text{III})} = 6.05 \pm 0.05$ in the resource-limited regime, to $m_\infty = 5.725 \pm 0.005$ otherwise. This mass reduction for abundant resource can be explained as follows: when volume-limited, cells keep dividing and vanishing after reaching confluence, hence they are in different elongation states. In contrast, resource-limited cells stop dividing

<footnote>$^4$With low initial cell count, as hysteresis exists at low resource.</footnote>
(they glide without oscillation) when they reach equilibrium, but they remain in an elongated state with enhanced mass below the division point.

The number of cells $v$, when sufficiently large, can be expressed in terms of the previously identified quantities in both regime $\frac{R}{|\Omega|} \ll \rho_*$ and $\frac{R}{|\Omega|} \gg \rho_*$:

$$v_* = \frac{\mu_\infty R}{m_* \rho_*} , \quad v_\infty = \frac{\mu_\infty |\Omega|}{m_\infty} , \quad (16)$$

with the specific relations between types I and II: $\rho_*(^I) = \rho_*(^II)$ and $m_*(^I) = r_*m_*(^II)$.

The equilibrium system with resource feedback effectively behaves as a re-tuned version of the original model without feedback. In the resource-limited regime $\frac{R}{|\Omega|} \ll \rho_*$, always reached in the large-volume limit $|\Omega| \to \infty$ at fixed $R$, the extensive behaviour of the reference model saturates at a total mass $\mu|\Omega| \sim \mu_\infty R/\rho_*$. During the transient growth, the system adjusts itself to a transition point on the critical manifold between the extensive (volume-limited) and non-extensive (resource-limited) phases of the parameter space. With feedback type II, all reference parameters are rescaled to $p/r_*$ by the common factor $r_* = 1 - \mu_\infty/\rho_*^{(I,II)} \approx 0.987$, while feedback type III only modifies the first parameter to $M_* = (1 + \mu_\infty/\rho_*^{(III)})M_c \approx 1.010 M_c$. Consequently, the hand-tuned reference parameters of Eq. (7) reside in the dense, volume-limited phase, within $\sim 1\%$ of the nearest point on the critical manifold.

## Summary and outlook

Strikingly, translating the Game of Life into a continuous and minimal model yields cell-like patterns with a nucleus and a shell, capable of division, motility, and death. We characterized their morphological and behavioural diversity through the lenses of bifurcations, symmetry breaking, shape instabilities, and critical dynamics. We refined its interpretation as a coarse-grained, non-standard reaction-diffusion system: its nonlinearity specifies homeostatic concentration ranges for morphogens, produced by a cascade of cell-sourced chemical species. Crucially, introducing a global conservation law and the associated coupling to a finite resource drives the system to self-organize at a point of phase transition between volume-limited and resource-limited dynamics. We now turn to perspectives and implications of this study.

Although the parameter space is reduced to only 6 + 1 degrees of freedom, finding a “living” region would be difficult without the history of successive rule refinements; the previous iteration was SmootherLife (CornusAmmonis, 2017). The extent of this region, especially for different scaling values $\lambda$, remains to be characterized. Continuous gliders and oscillators appear in $d = 1$, while cell division occurs in $d = 3$, but in very different regions of the parameter space, raising the question of how the spatial dimension $d$ shapes

cell morphology and parameter values. Such developments may uncover a more natural parametrization of the target function and lead to further parsimony.

Figure 3 suggests a dynamical phase transition, clarified by observing two regimes: abundant resource at fixed volume, and large volume at fixed resource. This approach is more general than the specific presence or absence of self-replicating patterns. A natural development consists in studying other signatures of criticality, such as power laws, as attempted in the discrete Game of Life (Alstrøm and Leão, 1994; Reia and Kinouchi, 2014).

As demonstrated with the global constraint of a conserved mass, certain types of environmental feedback act directly or indirectly as a homogeneous perturbation of a selected set of parameters. When introduced properly, this retroaction lets the system tune itself on the critical manifold (or close to it depending on the space and time limitations), effectively reducing the original parameter space by one degree of freedom. Including the tuning mechanism into the dynamics of existing models can accelerate the automated search for “living” regions close to phase transitions, as in (Papadopoulos et al., 2024).

Conservation laws also constitute a design principle for physically constrained, life-like systems. The local environmental retroaction of (Suzuki et al., 2023) constitutes a natural development. Its feedback mechanism via the resource abundance field $r(x, t)$ can be reformulated in terms of a diffusing available resource of concentration field $\bar{L} \propto r$, consumed or replenished when $L$ grows or decays: $\partial_t(L + \bar{L}) = \bar{D}\nabla^2\bar{L}$, which also satisfies global conservation Eq. (12) with $\bar{R} = \int \bar{L} dx$. The available resource is treated as a third fast-relaxing ($\bar{D} \to +\infty$) non-decaying auxiliary species, with uniform concentration $\bar{L}(t) = \bar{R}(t)/|\Omega|$ in our global approach. (Brauns et al., 2020) provide tools for analysing bifurcations and pattern formation in mass-conserving reaction-diffusion systems, while (Robinson et al., 2025) derive an effective Cahn-Hilliard-type coarse-grained description for large systems with locally conserved slow fields. These frameworks may inform the analysis and modelling of the non-standard nonlinearity $\Gamma$, which contains no spatial derivatives.

Finally, viewing the cGoL as a minimal model for hypothetical primordial life-forms raises the question of evolution. A biologically realistic evolutionary scenario would allow each cell to carry its own version of the parameters. It would open the possibility of individual mutations and their autonomous selection, propagation or rejection through the emergent behaviours of self-replication, motility, or death. The implementation of inter-individual diversity in the cGoL and its evolutionary application will be the topic of a future article.
Acknowledgements

We acknowledge the cross-fertilizing contributions from citizen science, generative arts and software engineering, which has fostered important progress, not only through models but also via efficient cross-platform implementations, interactive tools and empirical exploration techniques that help manage the complexity of these life-like systems and enable enthusiasts to explore their beauty. A.G. gratefully acknowledges Théotine Girardot, CornusAmmonis, Daniel M. Busiello, Jérémy Sourd and the anonymous reviewers for their valuable feedback and thoughtful suggestions, which played a significant role at various stages of development of this project. This work has been funded by the Max Planck Society.

References

Alstrøm, P. and Leão, J. (1994). Self-organized criticality in the "game of Life". *Physical Review E*, 49(4):R2507-R2508.

Brauns, F., Halatek, J., and Frey, E. (2020). Phase-Space Geometry of Mass-Conserving Reaction-Diffusion Dynamics. *Physical Review X*, 10(4):41036.

Calcaterra, C. and Boldt, A. (2023). Existence of life in Lenia. *International Journal of Dynamical Systems and Differential Equations*, 13(3):218-230.

Chan, B. W.-C. (2019). Lenia: Biology of Artificial Life. *Complex Systems*, 28(3):251-286.

Chan, B. W.-C. (2020). Lenia and Expanded Universe. In *The 2020 Conference on Artificial Life*, pages 221–229. MIT Press.

CornusAmmonis (2017). Smoother Life. https://www.shadertoy.com/view/XtVXzV.

Davis, Q. T. and Bongard, J. (2022). Step Size is a Consequential Parameter in Continuous Cellular Automata. In *The 2022 Conference on Artificial Life*. MIT Press.

Evans, K. M. (2001). Larger than Life: Digital Creatures in a Family of Two-Dimensional Cellular Automata. *Discrete Mathematics & Theoretical Computer Science*, DMTCS Proceedings vol. AA(Proceedings):2288.

Evans, K. M. (2003). Larger than Life: Threshold-range scaling of Life's coherent structures. *Physica D: Nonlinear Phenomena*, 183(1-2):45-67.

Gardner, M. (1970). Mathematical Games. *Scientific American*, 223(4):120-123.

Hutton, T., Munafo, R., Trevorrow, A., Rokicki, T., and Wills, D. (2015). Ready, a cross-platform implementation of various reaction-diffusion systems. https://github.com/GollyGang/ready.

Kawaguchi, T., Suzuki, R., Arita, T., and Chan, B. (2021). Introducing asymptotics to the state-updating rule in Lenia. In *The 2021 Conference on Artificial Life*. MIT Press.

Kojima, H. and Ikegami, T. (2023). Implementation of Lenia as a Reaction-Diffusion System. In *The 2023 Conference on Artificial Life*. MIT Press.

MacLennan, B. J. (1990). Continuous spatial automata.

Papadopoulos, V., Doat, G., Renard, A., and Hongler, C. (2024). Looking for Complexity at Phase Boundaries in Continuous Cellular Automata. In *Proceedings of the Genetic and Evolutionary Computation Conference Companion*, pages 179-182, Melbourne VIC Australia. ACM.

Pivato, M. (2007). *RealLife: The continuum limit of Larger than Life cellular automata*. *Theoretical Computer Science*, 372(1):46-68.

Plantec, E., Hamon, G., Etcheverry, M., Oudeyer, P.-Y., Moulin-Frier, C., and Chan, B. W.-C. (2023). Flow-lenia: Towards open-ended evolution in cellular automata through mass conservation and parameter localization. In *The 2023 Conference on Artificial Life*. MIT Press.

Rafler, S. (2011). Generalization of Conway's "Game of Life" to a continuous domain - SmoothLife. http://arxiv.org/abs/1111.1567.

Rafler, S. (2012). SmoothLife. https://sourceforge.net/projects/smoothlife/files/.

Reia, S. M. and Kinouchi, O. (2014). Conway's game of life is a near-critical metastable state in the multiverse of cellular automata. *Physical Review E*, 89(5):052123.

Rendell, P. (2002). Turing Universality of the Game of Life. In Adamatzky, A., editor, *Collision-Based Computing*, pages 513-539. Springer London, London.

Robinson, J. F., Machon, T., and Speck, T. (2025). Universal limiting behavior of reaction-diffusion systems with conservation laws. *Physical Review E*, 111(6):065417.

Suzuki, R., Asakura, K., and Arita, T. (2023). Lenia in a petri dish: Interactions between organisms and their environment in a Lenia with growth based on resource consumption. In *The 2023 Conference on Artificial Life*. MIT Press.