"Particle Lenia" 生命らしさを追求したライフゲーム
- NUM
- 3月21日
- 読了時間: 17分
更新日:4月13日
目次
はじめに
こんにちは NUMです。
今回はライフゲームの歴史とその進化の中で生まれた「Particle Lenia」という最新のシミュレーションモデルについてご紹介します。
僕は割と生物的なジェネラティブアート作品を作ることが多いです。自分のバックグラウンドも考えてみると、小さい頃から生き物は好きだったり、大学も生物系を出ていたりと根本に好きみたいな感情があるからそういう作品を作りがちなのかと思います。
最先端のライフゲーム「Particle Lenia」について仕組みを解説している日本語記事があまり見当たらなかったので、備忘録も兼ねて書いてみました。
この記事では「Particle Lenia」の数式とProcessingでの実装を対応させながら、なぜシンプルなルールから複雑なパターンが生まれるのかそしてParticle Leniaがどのような背景で誕生したのかを解説していきます。
ライフゲームの歴史
ライフゲームとは1970年にイギリスの数学者ジョン・ホートン・コンウェイ (John Horton Conway) が考案した数理モデルで、生命の四つの状態を使用してシンプルなルールでありながら予測不可能な複雑なパターンを生み出すことで大きな話題となったそうです。
この方の解説動画がわかりやすく面白いです。
基本ルール
各セルは生死の2状態を持ち、隣接するセルの状態に応じて次の状態が決まるというもので、わずかなルールの組み合わせで「生命」のような振る舞いが見られます。

継続的モデルへの展開
その後、多くの研究者やアーティストがライフゲームの魅力に着目し、離散的なセル・オートマトンを連続的なモデルに拡張する試みが始まりました。
元々のライフゲームは、格子状のセルが「生」か「死」の2値状態を取り、各セルが周囲のセルの状態に基づいて次の状態を決定するというシンプルなルールに基づいています。
しかし、実際の生物は徐々に「生」から「死」へ向かいます。 従来のライフゲームではこのような細かい変化や微妙なパターンの連続的な変化を表現するのが難しく、セルの状態が急激に変わるため、より自然な生命現象を再現するには限界がありました。
このような背景から、研究者たちはライフゲームのアイデアを「連続化」する試みを始めました。代表的な例が「SmoothLife」と「Lenia」です。
SmoothLife
SmoothLifeでは、セルの状態を0から1までの連続値として表現し、各セルの更新も連続的な数値として計算します。これにより、状態の変化が滑らかになり、より自然で有機的なパターンが生み出されるようになりました。
Lenia
Leniaはさらに進化したモデルで、連続空間と連続時間における更新ルールを導入しています。具体的には、各セルの状態更新に対して、周囲の状態の「畳み込み」を用いることで、局所的な影響をより正確に反映させ、生命体のような持続可能で変幻自在なパターンを生み出すことが可能になっています。
Particle Leniaとは
先ほど説明したLeniaは従来のグリッドベースモデルでしたが、Particle Leniaは粒子を個々の要素として扱うモデルです。
Particle Lenia Particleといえばシフマン先生の動画やNature of Codeでも出てきますが 基本的なオブジェクトの扱い方は同じで、粒子の振る舞いを計算する仕組みが少々複雑といったところです。
Particle Leniaには大きく2つの特徴があります。
引力と斥力 | 粒子同士は、引力(引き寄せる力) と 斥力(反発する力) のバランスによって動きます。どの粒子がどのくらい影響を及ぼすかは、「影響の範囲」と「影響の強さ」で決まります。
この影響の仕組みを調整することで、粒子は異なる動きを見せ、まるで生き物のような多様な模様や動きが生まれます。 |
エネルギー最小化 | エネルギー最小化とは、システム全体のエネルギーをできるだけ小さくしようとする原理で自然界でよく見られる振る舞いです。Particle Leniaでは次のような役割を果たします。
このプロセスにより、粒子たちは 「単なるランダムな動き」ではなく、「意味のある秩序を持つ動き」 を形成します。 |
数式
Particle Lenia は、連続空間上に存在する粒子同士が相互作用することで、生命現象に似た複雑なパターンが生じるモデルです。そのモデルを構成する数式について説明していきます。
カーネル関数 K(r)
各粒子が周囲に及ぼす影響は、まずカーネル関数によって表現されます。 カーネルは、粒子からの距離 r に応じた影響度を決める関数であり、一般的には次の形で定義されます。
カーネル関数で重要なのは以下式になります。
この式は正規分布の式で特定の範囲で影響を強くし、そこから離れると急激に影響を小さくする目的で使用されています。イメージは以下の通りです。

Lenia 場 U(x)
Lenia場は連続空間上の任意の点xにおける全粒子からの影響の合計です。
各粒子が点xに及ぼす影響(カーネル関数の値)を合算することで、その点の状態を求めます。
成長関数 G(u)
成長関数は粒子の成長効果や相互作用による引力を計算します。この関数により周囲の影響がある適正な範囲にある場合に、粒子が「引き寄せられる」効果が生まれます。
この数式の形、どこかで見たことありませんか?そうです正規分布の式です。 つまり、uの値がμ_Gに近い場合、Gの効果が大きくなり成長効果が大きく働くため、粒子同士が集まりやすくなります。 少しづつ数式が理解できるようになってきて楽しくなってきましたね!
反発場 R(x)
粒子同士があまりにも近づきすぎると、重なりを防ぐための反発効果が働きます。 具体的には「粒子との距離が 1 未満のときだけ急激にエネルギーが高まる」 という斥力の特性を持ちます。
反発場の数式で重要なのは以下式です。この数式は「xにいる粒子が、ほかの粒子pi とどれだけ近いか」を計算し、近すぎる場合に大きな反発エネルギーを発生させる役割があります。
特に 1 - ||x-p_i|| の数式がカギで、x、pi間の距離を1から引いている部分が近すぎる場合に大きな反発を生み出す部分です。本来、x、pi間の距離は1以上のになる場合がほとんどだからです。 ※今回はCrep/2やmax関数の二乗はプログラムをシンプルにする目的で省いています。
エネルギー場 E(x)
Particle Lenia は、自然界の自己組織化や生命現象を模倣するため
「エネルギーを下げるように粒子が動く」という仕組みを取り入れています。
全体のエネルギー E(x) は反発効果 R(x) と成長効果 G(U(x)) の差として以下数式で定義されます。 なぜ数式が引き算になるかというと、粒子が適度に集まると成長効果が強まり、結果的にエネルギーを下げる方向にさせるためです。
粒子の状態よって以下のようになると考えられます。
粒子が近づきすぎる
R(x) が増大し、エネルギーが上がる
粒子はその状態を避けるように離れようとする
粒子がほどよく集まる
G(U(x)) が大きくなり、エネルギーが下がる
粒子は適度な距離を保ちながら集団を作りやすくなる
粒子が遠く離れすぎる
G(U(x) は小さくなる(成長効果が得られない)
粒子はより引力が働く領域(他の粒子がそこそこ集まっている場所)へ移動しやすくなる
運動則:エネルギー勾配に沿った動き
最終的に、各粒子は自身が存在する位置でのエネルギー E(pi)を下げる方向へ動きます。 これは勾配降下法の考え方に基づいていて、運動則は以下のように表されます。
粒子 i の位置でのエネルギー勾配
各粒子はこの勾配の逆方向へ移動することで、システム全体のエネルギーを低下させ、より安定した状態を目指します。
実際のシミュレーションでは、オイラー法という数値積分法を用いて、離散的なタイムステップ dt ごとに次のように位置を更新します。
実装
final int N = 400; // 粒子数
final float mu_k = 4.0; // カーネル中心
final float sigma_k = 1.0; // カーネル幅
final float w_k = 0.022; // カーネル重み
final float mu_g = 0.6; // 成長関数中心
final float sigma_g = 0.15; // 成長関数幅
final float c_rep = 1.0; // 反発強度
final float dt = 0.1; // タイムステップ
Particle[] particles = new Particle[N]; // 粒子配列
void setup() {
size(500, 500);
pixelDensity(displayDensity());
frameRate(30);
// 初期は x,y ∈ [-2,2] の一様乱数で粒子を配置
for (int i = 0; i < N; i++) {
float x = random(-2, 2);
float y = random(-2, 2);
particles[i] = new Particle(new PVector(x, y));
}
}
void draw() {
background(0);
updateParticles(); // 各粒子の位置更新
for (Particle p : particles) {
p.show(5, 20);
}
}
/**
* 全粒子の相互作用に基づいてエネルギー勾配を計算し、
* Euler法を用いて各粒子の位置を更新
*/
void updateParticles() {
PVector[] newPositions = new PVector[N];
// 各粒子についてエネルギー勾配を計算
for (int i = 0; i < N; i++) {
Particle p_i = particles[i];
PVector gradR = calcRepulsiveGrad(p_i); // 反発効果の勾配 ∇R の計算
PVector gradG = calcGrowthGrad(p_i); // 成長効果の勾配 ∇G の計算
PVector gradE = PVector.sub(gradR, gradG); // エネルギー勾配 ∇E は ∇R - ∇G で求める
// Euler法による位置更新
newPositions[i] = PVector.sub(p_i.pos, PVector.mult(gradE, dt));
}
// 全粒子の位置を一括で更新
for (int i = 0; i < N; i++) {
particles[i].pos = newPositions[i];
}
}
/**
* カーネル関数 K(d) を計算
*
* @param 粒子間の距離
* @return カーネル関数の値
*/
float calcKernel(float d) {
return w_k * exp(-sq(d - mu_k) / sq(sigma_k));
}
/**
* カーネル関数 K(d) の微分 K'(d) を計算
*
* @param 粒子間の距離
* @return カーネル微分の値
*/
float calcKernelPrime(float d) {
return -2 * w_k * (d - mu_k) / sq(sigma_k) * exp(-sq(d - mu_k) / sq(sigma_k));
}
/**
* 指定した粒子 p_i に対する反発効果の勾配 ∇R を計算
*
* @param 計算対象の粒子
* @return 反発効果の勾配ベクトル
*/
PVector calcRepulsiveGrad(Particle p_i) {
PVector gradR = new PVector(0, 0);
// 他の全粒子との反発効果を計算
for (int j = 0; j < N; j++) {
Particle p_j = particles[j];
// 自身との相互作用は無視
if (p_i.equals(p_j)) continue;
// 粒子間の距離と方向を計算
PVector diff = PVector.sub(p_i.pos, p_j.pos);
float d = diff.mag();
// 非常に近い場合は計算をスキップ(数値安定性のため)
if (d < 1e-6) continue;
// 距離が 1 未満の場合に反発効果を計算
if (d < 1) {
PVector repTerm = diff.copy().normalize(); // 単位方向ベクトル
repTerm.mult(1 - d); // 距離に応じた反発強度
gradR.add(repTerm);
}
}
// c_rep を乗じ、負の方向に作用するので符号反転
gradR.mult(-c_rep);
return gradR;
}
/**
* 指定した粒子 p_i に対する成長効果の勾配 ∇G を計算
*
* @param 計算対象の粒子
* @return 成長効果の勾配ベクトル
*/
PVector calcGrowthGrad(Particle p_i) {
PVector gradU = new PVector(0, 0); // U の勾配
float U = 0; // 粒子の周囲の影響 U の総和
// 他の全粒子との相互作用を計算
for (int j = 0; j < N; j++) {
Particle p_j = particles[j];
// 自身との相互作用は無視
if (p_i.equals(p_j)) continue;
// 距離と方向を計算
PVector diff = PVector.sub(p_i.pos, p_j.pos);
float d = diff.mag();
if (d < 1e-6) continue;
// カーネル関数による影響を加算
float K_val = calcKernel(d);
U += K_val;
// カーネル微分による勾配を計算
float Kp = calcKernelPrime(d);
PVector termU = diff.copy().normalize(); // 単位方向ベクトル
termU.mult(Kp); // K'(d) * 単位ベクトル
gradU.add(termU);
}
float G = exp(-sq(U - mu_g) / sq(sigma_g)); // 成長関数 G(U) の計算
float dG_dU = -2 * (U - mu_g) / sq(sigma_g) * G; // 成長関数の微分
PVector gradG = gradU.copy().mult(dG_dU); // 連鎖律を適用した成長効果の勾配の計算
return gradG;
}
class Particle {
PVector pos; // 粒子の位置
/**
* Particle クラスのコンストラクタ
* @param 初期位置を表すPVector
*/
Particle(PVector pos) {
this.pos = pos;
}
/**
* 描画
*
* @param 描画時の粒子の直径
* @param シミュレーション座標から画面座標への変換スケール
*/
void show(float r, float scale) {
// シミュレーション座標から画面座標に変換
PVector screenPos = new PVector(pos.x * scale + width/2, pos.y * scale + height/2);
noStroke();
fill(255, 200);
circle(screenPos.x, screenPos.y, r);
fill(89, 183, 76, 100);
circle(screenPos.x, screenPos.y, r*2);
}
}
基本的にプログラムには詳細なコメントを記載しているので、重要な箇所を説明していきます。
void setup() {
size(500, 500);
pixelDensity(displayDensity());
frameRate(30);
// 初期は x,y ∈ [-2,2] の一様乱数で粒子を配置
for (int i = 0; i < N; i++) {
float x = random(-2, 2);
float y = random(-2, 2);
particles[i] = new Particle(new PVector(x, y));
}
}
void draw() {
background(0);
updateParticles(); // 各粒子の位置更新
for (Particle p : particles) {
p.show(5, 20);
}
}
プログラムの単純化すると3つの処理になります。
Particleをインスタンス化して初期配置を決定
Particleの位置更新
Particleの描画
class Particle {
PVector pos; // 粒子の位置
/**
* Particle クラスのコンストラクタ
* @param 初期位置を表すPVector
*/
Particle(PVector pos) {
this.pos = pos;
}
/**
* 描画
*
* @param 描画時の粒子の直径
* @param シミュレーション座標から画面座標への変換スケール
*/
void show(float r, float scale) {
// シミュレーション座標から画面座標に変換
PVector screenPos = new PVector(pos.x * scale + width/2, pos.y * scale + height/2);
// 微生物に近づけるために緑と白の円を描画
noStroke();
fill(255, 200);
circle(screenPos.x, screenPos.y, r);
fill(89, 183, 76, 100);
circle(screenPos.x, screenPos.y, r*2);
}
}
この記事ではParticle Leniaの動作原理に焦点を当てるため、Particleクラスの構造もシンプルに位置メンバと描画メソッドのみとしています。
/**
* 全粒子の相互作用に基づいてエネルギー勾配を計算し、
* Euler法を用いて各粒子の位置を更新
*/
void updateParticles() {
PVector[] newPositions = new PVector[N];
// 各粒子についてエネルギー勾配を計算
for (int i = 0; i < N; i++) {
Particle p_i = particles[i];
PVector gradR = calcRepulsiveGrad(p_i); // 反発効果の勾配 ∇R の計算
PVector gradG = calcGrowthGrad(p_i); // 成長効果の勾配 ∇G の計算
PVector gradE = PVector.sub(gradR, gradG); // エネルギー勾配 ∇E は ∇R - ∇G で求める
// Euler法による位置更新
newPositions[i] = PVector.sub(p_i.pos, PVector.mult(gradE, dt));
}
// 全粒子の位置を一括で更新
for (int i = 0; i < N; i++) {
particles[i].pos = newPositions[i];
}
}
updateParticles関数はParticle Leniaの肝となる関数になります。 この関数は粒子一つずつの位置を計算し、次の位置を決定します。 この後出てきますが反発効果(calcRepulsiveGrad関数)、成長効果(calcGrowthGrad関数)は自身と自身以外の粒子の近さを考慮することで計算ができます。 (現段階では全体理解のため詳細は省きます。) 先ほど計算した反発効果、成長効果の値を使用してエネルギー勾配を求めます。
PVector gradR = calcRepulsiveGrad(p_i); // 反発効果の勾配 ∇R の計算
PVector gradG = calcGrowthGrad(p_i); // 成長効果の勾配 ∇G の計算
PVector gradE = PVector.sub(gradR, gradG); // エネルギー勾配 ∇E は∇R-∇Gで求める
オイラー法を使用して次の位置を計算します。 オイラー法は現在の状態で計算した勾配(変化率)をそのまま使って次の状態を求める方法です。
// Euler法による位置更新
newPositions[i] = PVector.sub(p_i.pos, PVector.mult(gradE, dt));
/**
* 指定した粒子 p_i に対する反発効果の勾配 ∇R を計算
*
* @param 計算対象の粒子
* @return 反発効果の勾配ベクトル
*/
PVector calcRepulsiveGrad(Particle p_i) {
PVector gradR = new PVector(0, 0);
// 他の全粒子との反発効果を計算
for (int j = 0; j < N; j++) {
Particle p_j = particles[j];
// 自身との相互作用は無視
if (p_i.equals(p_j)) continue;
// 粒子間の距離と方向を計算
PVector diff = PVector.sub(p_i.pos, p_j.pos);
float d = diff.mag();
// 非常に近い場合は計算をスキップ(数値安定性のため)
if (d < 1e-6) continue;
// 距離が 1 未満の場合に反発効果を計算
if (d < 1) {
PVector repTerm = diff.copy().normalize(); // 単位方向ベクトル
repTerm.mult(1 - d); // 距離に応じた反発強度
gradR.add(repTerm);
}
}
// c_rep を乗じ、負の方向に作用するので符号反転
gradR.mult(-c_rep);
return gradR;
}
続いて反発効果を計算するcalcRepulsiveGrad関数になります。 この関数はcalcGrowthGrad関数と構造が似通っていて、計算対象の粒子とそれ以外の全粒子の距離を比較させて相互作用の具合を計算するので、全粒子をループさせています。 以下は共通して存在するスキップ処理です。
// 自身との相互作用は無視
if (p_i.equals(p_j)) continue;
// 粒子間の距離と方向を計算
PVector diff = PVector.sub(p_i.pos, p_j.pos);
float d = diff.mag();
// 非常に近い場合は計算をスキップ(数値安定性のため)
if (d < 1e-6) continue;
// 距離が 1 未満の場合に反発効果を計算
if (d < 1) {
PVector repTerm = diff.copy().normalize(); // 単位方向ベクトル
repTerm.mult(1 - d); // 距離に応じた反発強度
gradR.add(repTerm);
}
粒子間の距離d が 1 未満のときに反発力を計算しています。粒子が近すぎると「押し返す力」を働かせたいためです。
反発を働かせる方向を単位方向ベクトルを計算することで取得し、距離に応じた反発強度にスケールしています。(近いほど強い反発を与え、遠いほど弱くするための係数)
そして、全粒子の相互作用を加味した反発を計算するためgradR.add(repTerm);を実施します。
Eの式はエネルギー最小化が目的となるため最終的な更新式でマイナス符号を適用します。
// c_rep を乗じ、負の方向に作用するので符号反転
gradR.mult(-c_rep);
数式のセクションでこのマイナスについては言及されていなかった理由は力の計算原理に焦点を置いていたためです。ただ、プログラムでは 実際にどの方向に動かすかをはっきりさせる必要がある ので、gradR.mult(-c_rep);のようにマイナスを付与して方向を反転させています。
/**
* 指定した粒子 p_i に対する成長効果の勾配 ∇G を計算
*
* @param 計算対象の粒子
* @return 成長効果の勾配ベクトル
*/
PVector calcGrowthGrad(Particle p_i) {
PVector gradU = new PVector(0, 0); // U の勾配
float U = 0; // 粒子の周囲の影響 U の総和
// 他の全粒子との相互作用を計算
for (int j = 0; j < N; j++) {
Particle p_j = particles[j];
// 自身との相互作用は無視
if (p_i.equals(p_j)) continue;
// 距離と方向を計算
PVector diff = PVector.sub(p_i.pos, p_j.pos);
float d = diff.mag();
if (d < 1e-6) continue;
// カーネル関数による影響を加算
float K_val = calcKernel(d);
U += K_val;
// カーネル微分による勾配を計算
float Kp = calcKernelPrime(d);
PVector termU = diff.copy().normalize(); // 単位方向ベクトル
termU.mult(Kp); // K'(d) * 単位ベクトル
gradU.add(termU);
}
float G = exp(-sq(U - mu_g) / sq(sigma_g)); // 成長関数 G(U) の計算
float dG_dU = -2 * (U - mu_g) / sq(sigma_g) * G; // 成長関数の微分
PVector gradG = gradU.copy().mult(dG_dU); // 成長効果の勾配の計算
return gradG;
}
続いて成長効果(粒子を引き寄せる引力)を計算するcalcRepulsiveGrad関数になります。
Gの式は内部にUの式を内包しているので連鎖律を考慮した計算式にする必要があります。
連鎖律の適用により内側の関数Uの勾配と外側の関数Gの微分を掛け合わせることで成長効果の計算を実現しています。
PVector gradU = new PVector(0, 0); // U の勾配
float U = 0; // 粒子の周囲の影響 U の総和
// 他の全粒子との相互作用を計算
for (int j = 0; j < N; j++) {
Particle p_j = particles[j];
// 自身との相互作用は無視
if (p_i.equals(p_j)) continue;
// 距離と方向を計算
PVector diff = PVector.sub(p_i.pos, p_j.pos);
float d = diff.mag();
if (d < 1e-6) continue;
// カーネル関数による影響を加算: U = Σ K(d)
float K_val = calcKernel(d);
U += K_val;
// カーネル微分による勾配を計算: 単位ベクトル * K'(d)
float Kp = calcKernelPrime(d);
PVector termU = diff.copy().normalize(); // ここで単位方向ベクトルを取得
termU.mult(Kp); // K'(d) を掛ける
gradU.add(termU); // これを全粒子分足し合わせる
}
まずは成長関数の微分で使用するためのU(x)と∇Uを前段のforループで自身以外の粒子を考慮して計算しています。可読性を重視してカーネル関数とカーネル関数の微分計算は以下に外出ししています。勾配の計算はcalcRepulsiveGrad関数と仕組みは同じなので割愛します。
/**
* カーネル関数 K(d) を計算
*
* @param 粒子間の距離
* @return カーネル関数の値
*/
float calcKernel(float d) {
return w_k * exp(-sq(d - mu_k) / sq(sigma_k));
}
/**
* カーネル関数 K(d) の微分 K'(d) を計算
*
* @param 粒子間の距離
* @return カーネル微分の値
*/
float calcKernelPrime(float d) {
return -2 * w_k * (d - mu_k) / sq(sigma_k) * exp(-sq(d - mu_k) / sq(sigma_k));
}
次は成長関数の微分になります。先ほど計算したU(x)はここで使用されていますね。
とはいっても特に難しいことはなく、この式は単純にプログラムに落とし込むことが可能なので問題ないと思います。事前にGを計算し、dG_dUを求めているだけです。 最後のgradU.copy().mult(dG_dU)は∇G(U(x))の全体式であった∇Uを掛けている部分に対応する箇所になります。
float G = exp(-sq(U - mu_g) / sq(sigma_g)); // 成長関数 G(U) の計算
float dG_dU = -2 * (U - mu_g) / sq(sigma_g) * G; // 成長関数の微分
PVector gradG = gradU.copy().mult(dG_dU); // 成長効果の勾配の計算
まとめ
今回はParticle Leniaの数式原理と実装を解説してみました。
ライフゲームには長い歴史があり、数々の研究者が生命らしさの記述を試みた結果、Particle Leniaというモデルが生まれたんだと実感しました。ALIFEにますますハマりそうです、、
また今回の記事は原理を理解するのにかなり苦労しましたが、以前書いたチューリングパターンの記事や今までの知識が使える場面が多く、実感はないですが力はついているんだなと感じました。
今までの傾向から生命現象あるあるを言うと、
「勾配」というワードが出てきがち。

この作品はParticle Leniaを使用して、分化前の細胞をイメージして制作しました。
分化とは細胞分裂を繰り返していくと、特定の機能を持つように分化していきます。 例えば脳、目、口とかですね。
分化前でありつつも様々な可能性を秘めていると言う意味を込めて、いつも使用する赤、青、黄の細胞を沢山配置し、力を発散する寸前という表現をしました。 要は「何にでもなれる」状態をビジュアル化したような感じです。 最後まで読んでいただきありがとうございました!
参考
Comments