| Authors | Affiliation | Published |
|---|---|---|
| Alexander Mordvintsev | December 23, 2022 | |
| Eyvind Niklasson | ||
| Ettore Randazzo | Reproduce in Notebook |
Patterns and behaviors you see in the video above were produced by the artificial life system that we called âParticle Leniaâ. Itâs inspired by Lenia Mathematical Lifeforms, which in turn is a continuous generalization of John Conway's Game of Life. In this article we are going to describe the motivation, implementation and behavior of Particle Lenia.
Before diving into the details of the proposed model, we would like to speculate a little about why to create such models. Physics, chemistry and biology study the organization of matter at various scales, from subatomic to cosmic. The behavior of each scale is determined by the laws that govern behavior of the underlying scale, but the nature of relations between scales is full of surprises. For example letâs consider a simple idea: matter consists of tiny, constantly moving particles, which attract each other at close distances, but repel when they are too close. Lennard-Jones potential is one of possible mathematical formalizations of the above idea. In spite of the apparent simplicity of microscopic dynamics, at large scales these rules give rise to a plethora of complex phenomena, such as aggregate states of matter, thermodynamics, acoustics, hydrodynamics and turbulence. In scientific and engineering practice we can often abstract microscopic details behind macroscopic models, such as diffusion or Navier-Stokes equations, that operate at the scale of interest. But in some domains, like material science or biochemistry, we have to focus on the relationship between many orders of magnitude in scale. How do changes in molecular structure reflect in macroscopic material properties? How do genomic sequences translate into morphological appearance? In the field of biology these ties between scales are especially important and obscure at the same time. Thatâs why development of practical methods of capturing and analyzing such relations are of a great need. In silico multi-agent models provide a computational playground for creation and improvement of such methods. We may use such models for efficient design of local microscopic rules that lead to a desired global behavior. Thus, creation of simulated Artificial Life serves as a stepping stone in deep understanding of life in the real world.
Game of Life cellular automata captivated the attention of a few generations of researchers and hobbyists with a diversity of complex behaviors stemming from a very simple set of local rules. A few continuous generalizations were proposed, SmoothLife and Lenia are most popular among them. They describe evolution of one or a few scalar fields with a system of PDEs. Such models can support a great diversity of interesting and beautiful life-like behaviors. Yet we argue that PDE-based models have a couple of drawbacks, which motivated the creation of the particle based analogy.
First, SmoothLife and Lenia âcreaturesâ are compact persistent patterns that are defined by both spatial configuration of the field (âmorphologyâ), and parameters of the PDE system that governs its evolution (âlaws of physicsâ). Thus, two creatures can coexist in the same universe only when they are supported by the same set of parameters. Otherwise, some mechanism is required that would modify local simulation rules depending on the spatial field configuration.
Second, in their original formulation neither SmoothLife nor Lenia have any conservation laws, which leads to often observed explosive growth or extinction behaviors.
Recently proposed Flow Lenia model tries to work around these issues by formulating the PDE system around a flow field that moves the distribution of âmatterâ in space without creation of a new mass. Authors also propose a method of combining multiple update rules in the same space. Particle Lenia can be seen as another perspective on the Flow Lenia model: we represent moving matter as a population of particles, rather than a scalar field (see Lagrangian flow vs. Eulerian).
We think that particle based representation has a number of attractive properties:
Mass conservation simply follows from the constant number of particles. Grid-based methods need to take special measures to accurately redistribute mass between grid cells.
Combining multiple different update rules is possible by attaching them to individual particles (note that we mostly focus on single particle-type models in this article, as we find their behavior sufficiently interesting).
Particle based simulation is easy to implement and scales well when switching from 2D to 3D space.
There are quite a few known particle-based artificial life simulations, including the classic Boids, Primordial Particle Systems and many others. We think that the main contribution of this work is to provide a link between the popular Lenia model and particle systems. The proposed artificial life model demonstrates a large diversity of interesting and even beautiful behaviors, even in the case when all particles share the same parameters. We also examine the model dynamics through the lens of a greedy local energy minimization process, and compare the behavior of the model to that of running with an equivalent global energy minimization routine. We find that the difference in behavior between these paradigms is non-intuitive and may have significance for the design of future systems of optimized interacting agents.
Many classic physics models describe particles that interact through fields they induce. Canonical examples are point charges and electrostatic field, point masses and gravity, or aforementioned Lennard-Jones particles. Particle positions define the configuration of the potential field, and the field affects the motion of particles, which accelerate towards the gradient of the potential field. Particle Lenia follows the same principle. In this section we are going to develop and implement Particle Lenia step by step, while contrasting it against the original Lenia model.
Lenia is a PDE system that describes the evolution of the scalar field $\mathbf{A}$:
$$\mathbf{A}^{t+\Delta t} = \left[ \mathbf{A}^{t} + \Delta t \;G \big(\mathbf{K}*\mathbf{A}^t\big) \right]_0^1$$
The behavior of the system is defined by the initial state of the field $\mathbf{A}^0$, the scalar growth mapping $G$ and the convolution kernel $\mathbf{K}$. The field $\mathbf{U}^t=\mathbf{K}*\mathbf{A}^t$ is referred to as the potential distribution in the original Lenia paper, and $\mathbf{G}^t=G \big(\mathbf{U}^t\big)$ is called the growth distribution. In this work we are going to use a different naming convention to avoid confusion in the discussion about energy-based interpretation that will follow later. We are going to call $\mathbf{G}^t$ â growth field and $\mathbf{U}^t$ â Lenia field.
Lenia field. Particle Lenia is an ODE system in which the scalar field $\mathbf{A}$ is replaced by a population of particles, characterized by their location vectors $\mathbf{p}_i^t$. We can still compute the value of the continuous field $\mathbf{U}^t$ at any point $\mathbf{x}$ by treating each particle as a unit impulse and adding contributions of individual particles:
$$\mathbf{U}^t(\mathbf{x}) = \sum_i \mathbf{K}(\mathbf{x}-\mathbf{p}_i^t) = \sum_i K\big(\Vert\mathbf{x}-\mathbf{p}_i^t\Vert\big)$$
$$K(r) = w_K\exp\big({-(r-\mu_K)^2 \big/ \sigma_K^2}\big)$$
As in Lenia, we choose $\mathbf{K}$ to be radially symmetrical, depending only on the distance from the particle. This symmetry makes the model rotationally equivariant. In this work we use the simplest Lenia kernel, that consists of a single circular (2d) or spherical (3d) shell. The coefficient $w_K$ is selected in such a way that the integral of $\mathbf{K}$ over the whole space equals one. Growth field. Knowing $\mathbf{U}^t$ we can compute $\mathbf{G}^t=G \big(\mathbf{U}^t\big)$. Like in Lenia, $G$ is a single gaussian peak that "selects" soft range of $\mathbf{U}^t$ values where growth is happening, or in our case, particles are attracted. Originally $G$ had values in $[-1,1]$ range, but in our case, as we show later, only its derivative matters, so additional scaling isn't needed:
$$G(u) = \exp\big({-(u-\mu_G)^2 \big/ \sigma_G^2}\big)$$
Repulsion field. Lenia explicitly bounds values of $\mathbf{A}$ to stay in $[0, 1]$ range. We decided to implement something analogous to this. The lower boundary is trivial because there can't be less than zero particles in any area. We added a repulsion force between particles to limit the maximum particle concentration. The magnitude of force that repels two particles $i \neq j$ is $f_{ij} = f_{ji}= c_\text{rep}\max(1-\Vert\mathbf{p}_i-\mathbf{p}_j\Vert, 0)$, where $c_\text{rep}$ controls the repulsion strength. We can also represent this force with a repulsion potential field $\mathbf{R}^t$:
$$\mathbf{R}^t(\mathbf{x}) = \frac{c_\text{rep}}{2}\sum_{i: \mathbf{p}_i \neq \mathbf{x}} \max\big(1-\Vert \mathbf{x} - \mathbf{p}_i\Vert,\; 0\big)^2$$
The gradient of this field equals to the sum of repulsion forces acting on a particular point in space. Motion law. Particles in our model move against the local gradient of the energy field $\mathbf{E}=\mathbf{R}-\mathbf{G}$. In the language of ODEs it can be simply formulated as:
$$ \frac{d\mathbf{p}_i}{dt} = -â\mathbf{E}(\mathbf{p}_i) = -\left[ \frac{â\mathbf{E}(\mathbf{p}_i)}{â\mathbf{p}_i} \right]^\top $$
Particles are trying to minimize repulsion energy $\mathbf{R}$, while seeking for areas with high values for $\mathbf{G}$. Note that they don't have momentum and directly follow in the direction of the applied force. We are trying to get the dynamics similar to the behavior of microscopic organisms in water, where kinetic energy dissipates almost immediately if force not applied.
Let's express the formulas above with JAX code:
Params = namedtuple('Params', 'mu_k sigma_k w_k mu_g sigma_g c_rep')
Fields = namedtuple('Fields', 'U G R E')
def peak_f(x, mu, sigma):
return jp.exp(-((x-mu)/sigma)**2)
def fields_f(p: Params, points, x):
r = jp.sqrt(jp.square(x-points).sum(-1).clip(1e-10))
U = peak_f(r, p.mu_k, p.sigma_k).sum()*p.w_k
G = peak_f(U, p.mu_g, p.sigma_g)
R = p.c_rep/2 * ((1.0-r).clip(0.0)**2).sum()
return Fields(U, G, R, E=R-G)
def motion_f(params, points):
grad_E = jax.grad(lambda x : fields_f(params, points, x).E)
return -jax.vmap(grad_E)(points)
That's it. fields_f function, given the system parameters and positions of particles, computes values of all fields described earlier at the point x. motion_f combines grad and vmap to compute velocities of all points at a given system state.
Now we can simply pick some parameter values, sample random points and apply and ODE integrator (even Euler would work):
params = Params(mu_k=4.0, sigma_k=1.0, w_k=0.022, mu_g=0.6, sigma_g=0.15, c_rep=1.0)
key = jax.random.PRNGKey(20)
points0 = (jax.random.uniform(key, [200, 2])-0.5)*12.0
dt = 0.1
def odeint_euler(f, params, x0, dt, n):
def step_f(x, _):
x = x+dt*f(params, x)
return x, x
return jax.lax.scan(step_f, x0, None, n)[1]
rotor_story = odeint_euler(motion_f, params, points0, dt, 10000)
animate_lenia(params, rotor_story, name='rotor.mp4')
Left plot shows lenia field $\mathbf{U}$ in blue, with yellow overlay $\mathbf{G}$. Right visualization shows the evolution of the "energy" field $\mathbf{E}$. White color corresponds to the empty space energy level equal to $G(0)$. Red areas show higher energy levels produced by the repulsion field $\mathbf{R}$. In green areas the growth field $\mathbf{G}$ outweighs repulsion. The black dash on the colorbar corresponds to the average value of $\mathbf{E}$ at current particle locations.
We see these simple rules can produce complex and interesting behaviors. The system undergoes a sequence of transitions through complex formations, some of which look almost stable and ends up in a spinning attractor. Different random seeds produce different scenarios and sometimes even lead to a different final state.
The total system energy that we express as $E_\text{total}=\sum_i \mathbf{E}(\mathbf{p}_i)$ provides an interesting perspective on the system dynamics:
def point_fields_f(params, points):
return jax.vmap(partial(fields_f, params, points))(points)
def total_energy_f(params, points):
return point_fields_f(params, points).E.sum()
As we see, the total system energy is mostly decreasing through a series of phase transitions between different formations. Nevertheless there are intervals where the energy is increasing. Why is this happening? Shouldn't the system that follows the energy gradient, eventually find a stable local minima? Turns out that there is one subtle detail which makes the dynamics much more interesting. Each particle greedily minimizes its own energy without considering its influence on other particles. We can propose a different motion rule, where all particles move in the direction that minimizes the total energy:
$$ \frac{d\mathbf{p}_i}{dt} = - \left[ \frac{âE_\text{total}}{â\mathbf{p}_i} \right]^\top = - \left[ \frac{â \sum_j \mathbf{E}(\mathbf{p}_j)}{â\mathbf{p}_i} \right]^\top$$
We can easily express this intent using automatic differentiation:
def total_motion_f(params, points):
return -jax.grad(partial(total_energy_f, params))(points)
In the following video we run the global energy descent for the first 2000 steps, then switch to the local descent. We can see that the global optimization indeed drives the system into a stable low-energy state. Once we switch to the local rule, the system escapes this state and enters a dynamic attractor similar to the one from the previous simulation.
If we compare energy plots from two simulations, we see that the global descent finds a much lower energy state that can't be sustained by the local update rule. Nevertheless we focus on the local descent rule because we think that it produces much more interesting and diverse structures and behaviors. Even more surprisingly, it was recently shown (EigenGame) that in some cases multi-agent local optimization opens new perspectives on solving well known computational problems, such as PCA.
Energy evolution analysis opens ways to quickly compare and spot patterns and irregularities in large ensembles of Particle Lenia simulations. For example, the figure below shows energy histories of 100 simulations that only differ by random seed used to sample starting point positions. We also show final configurations of four selected runs.
Original Lenia is a PDE system, which means that in order to simulate it on a computer we have to discretize the space somehow. It was shown that Lenia creatures can exist on regular and even non-regular grids of various resolutions. Of course, discretization artifacts can't be completely neglected, for example some behaviors can't exist if the grid resolution is too low, while other behaviors may exploit the raster grid structure and exist at low resolutions only.
In Particle Lenia, the number of particles can be seen as an analog to the grid resolution. In this section we try to see if it's possible to reproduce a specific pattern with a different number of particles. Consider a 200-particle "rotor" pattern we discovered in the previous section (third from left on the figure above). Let's apply following transformations to it:
Scale point coordinates and extent of the kernel $\mathbf{K}$ by two: $\mathbf{p}'_i = 2\mathbf{p}_i, \mu_K'=2\mu_K, \sigma_K'=2\sigma_K$. It's easy to see that this transformation doesn't change the values of $\mathbf{U}$ (and thus $\mathbf{G}=G(\mathbf{U})$) perceived by particles, but they are distributed much sparser now, so repulsion forces have weakened or vanished.
Replace each particle with four new ones to restore density and repulsion force. The previous step increased the area covered by the formation four times. Now we replace each particle with four new ones, distributed around the old particle position in a rotated cross fashion (similar to the rotated grid super-sampling)
Fix the range of $\mathbf{U}$. The previous step quadrupled the number of particles, which in turn made values of $\mathbf{U}$ four times larger on average. To compensate, we may either set $w_K' = w_K/4$, or adjust $\mu_G'=4\mu_G, \sigma_G'=4\sigma_G$. We prefer to use the former method.
The above procedure is implemented in the following code. Let's try to upscale the third pattern from the previous diagram:
def upscale_lenia(p: Params, points):
params2 = p._replace(mu_k=p.mu_k*2, sigma_k=p.sigma_k*2, w_k=p.w_k/4)
offsets = jp.float32([[3,1], [-1,3], [-3,-1], [1,-3]])/6
points2 = (points[:,None,:]*2 + offsets).reshape(-1, 2)
return params2, points2
params2, points2 = upscale_lenia(params, stories[52,-1])
story2 = odeint_euler(motion_f, params2, points2, dt*2, 10000)
animate_lenia(params2, story2, rate=20, w=600, show_UG=False)
The proposed method allowed us to reproduce the 200-particle "rotor" behavior with 800 particles. Note that in general we are not guaranteed to succeed, as some patterns may depend on the granularity of the "matter" to function properly. We also observed that the "rotor" is quite sensitive to small perturbations and can easily fall into a lower energy attractor.
You may notice that the set of parameters we use is a bit redundant. In particular, multiplying $w_K$ and dividing $\mu_G$ and $\sigma_G$ by the same coefficient doesn't have an impact on the system dynamics. To keep $\mu_G$ and $\sigma_G$ in the same range across different scales, we use numerical integration to normalize the kernel $\mathbf{K}$:
def calc_K_weight(mu, sigma, dim_n):
r = jp.linspace(max(mu-sigma*4, 0.0), mu+sigma*4, 51)
y = peak_f(r, mu, sigma) * r**(dim_n-1)
s = jp.trapz(y, r) * {2:2,3:4}[dim_n] * jp.pi
return 1.0/s
def create_params(m_k, s_k, m_g, s_g, rep, dim_n):
w_k = calc_K_weight(m_k, s_k, dim_n)
return Params(m_k, s_k, w_k, m_g, s_g, rep)
One of the benefits of particle-based approaches is the ease of scaling to higher dimensions. The Lenia code we wrote above can simulate behavior of 3D particles without modification. By the way, 3D rendering is implemented in pure JAX, please refer to this article for details.
def sample_ball(key, n, R=8.4, dim_n=3):
k1, k2 = jax.random.split(key)
pos = normalize(jax.random.normal(k1, [n, dim_n]))
return pos * jax.random.uniform(k2, [n, 1])**(1/dim_n)*R
params3d = create_params(4.3, 1.1, 0.23, 0.06, 1.0, dim_n=3)
key = jax.random.PRNGKey(13)
points = sample_ball(key, 256)
story3d = odeint_euler(motion_f, params3d, points, 0.2, 15000)
with VideoWriter() as vid:
for i, p in enumerate(story3d[::20]):
vid(render_lenia3d(params3d, p,
camera_pos=jp.array([jp.cos(i/100), 0.5, jp.sin(i/100)])*13.))
This video highlights another aspect of Particle Lenia. Particles are initially uniformly scattered in space. Then a cluster with organized structure starts to grow. The shape and behavior of this cluster dramatically changes as it absorbs more particles. Eventually particles assemble into a mushroom-shaped glider that slowly drifts away. It seems that this shape is an attractor for a given set of Lenia parameters and number of points, as running simulation with different initial conditions produces the similar gliders, flying in different directions.
While working on this article we experimented with various ways of visualizing particles and fields surrounding them. We also thought that it might be interesting to create a kind of audio representation of the particle system state. The process of turning data into sounds is often called sonification. We tried a simple method when every particle is given a voice, which has a frequency depending on the particle potential energy and volume is proportional to the particle current speed. The video below shows an example sonification of the simulation from the beginning of this article.
The text above contains a number of references to other works, which we would like to summarize and categorize in this section.
Lenia and Flow Lenia Particle Lenia was directly inspired by Flow Lenia (arxiv, colab, videos) and is designed with the same goals in mind (mass conservation and multiple species in the same world). Our model description closely follows terminology and notation used in original Lenia paper (arxiv, video). We would also like to mention the Lenia Expanded Universe (arxiv, video) and Sensorimotor Lenia articles, because they provide a clear path for future development of our work.
Energy-based models The energy-based formulation of Particle Lenia is inspired by a number of classic models from physics and chemistry. Lennard Jones potential is a classic example of a system that combines attractive and repulsive potential fields. Energy terms of the Embedded Atom Model have even greater similarity to the Particle Lenia energy. The biggest difference of our model is the lack of momentum; the assumption that particles are effectively moving through a viscous substance that quickly dissipates kinetic energy. We also think it is worth mentioning JAX MD that implements a number of molecular models using JAX differentiable programing framework.
Particle-based artificial life Artificial Life literature is full of examples of simple muliagent models that have sophisticated behaviors. This phenomena is often refered to as emergence, self-organization or swarm intelligence. Boids model of flocking is a classic example from that category. There are plenty other interesting examples, such as Primordial Particle Systems and models based based on non-reciprocal forces between different particle types. Other models try to mimic pheromone signaling through the environment (e.g. Physarum model).
Multi-agent optimization Particle Lenia model can be interpreted as a collective of agents that try to minimize their individual objective functions. This gives rise to behaviors that have dramatically different characteristics from global energy minimization. Recent work on EigenGame provides an interesting example when local multi-agent optimization has a number of desirable traits in comparison to global optimization. These traits enable new ways of solving classic computational problems.
In this article we described Particle Lenia, a new particle system ODE-based artificial life model that is inspired by the Lenia PDE system.
Interactive demo and collecting creatures. The introduction video shows a few creatures and behaviors we found by varying Particle Lenia parameters. These creatures were discovered using the interactive browser-based demo. We invite readers to participate in finding new creatures and sharing them on social media. Preview images created by the demo contain all information needed to bring them into life again, simply drag and drop the image file into the demo window!
Future work The obvious bit that is missing from this article are multi-creature/multi-particle types experiments. Although the model simulation code can be trivially extended to support per-particle parameters, we decided to leave this exploration for the future. Another promising future work direction is goal driven search for patterns and rules capable of forming specific patterns or performing specific tasks, in the spirit of Sensorimotor Lenia. We believe that these principles of particle swarm design can later be applied to real world self-organizing systems.
Reproducibility All experiments and diagrams in this article were produced using this colab notebook. The introduction video demonstrates a number of creatures that were discovered, and can be explored in the supplementary interactive demo
Acknowledgements We would like to thank Bert Chan and INRIA Flowers team for inspiring this work and providing early feedback.
Citing Particle Lenia For attribution in academic contexts, please cite the arxiv version of this article (to be released Q1 2023)