パターンの探索 - 数理生物学演習

パラメータや初期条件を変え，反応拡散モデルにどんなパターンが現れるかを探る．

このページを単独で実行できるよう，前節の関数をまとめて再掲する．

import matplotlib.pyplot as plt
import numpy as np


def laplacian(f, dh):
    """周期境界条件でのラプラシアン（ベクトル化）."""
    return (
        np.roll(f, 1, axis=0)
        + np.roll(f, -1, axis=0)
        + np.roll(f, 1, axis=1)
        + np.roll(f, -1, axis=1)
        - 4 * f
    ) / dh**2


def update(u, v, D_u, D_v, mu_u, mu_v, k_1, k_2, k_3, dh, dt):
    """反応拡散モデル（ギーラー-マインハルト系）の1ステップ更新（周期境界条件）."""
    u_new = u + (D_u * laplacian(u, dh) - mu_u * u + k_1 * u**2 / v + k_2) * dt
    v_new = v + (D_v * laplacian(v, dh) - mu_v * v + k_3 * u**2) * dt
    return u_new, v_new

パラメータを変えて何度も実行するので，シミュレーションを1つの関数にまとめておく．拡散係数 D_u 以外は基本値を既定値にしておき，変えたいものだけを指定できるようにする．

dt = 0.02  # 時間の刻み幅
dh = 1.0  # 空間の刻み幅
x = np.arange(-50, 50, dh)
y = np.arange(-50, 50, dh)


def simulate(
    D_u, D_v=10.0, mu_u=1.0, mu_v=1.0, k_1=1.0, k_2=0.05, k_3=1.0, t_end=1000.0
):
    """指定したパラメータで反応拡散シミュレーションを実行し，最終状態の u を返す."""
    rng = np.random.default_rng()
    u = 1.0 + rng.uniform(-0.001, 0.001, size=(len(y), len(x)))
    v = 1.0 + rng.uniform(-0.001, 0.001, size=(len(y), len(x)))
    for _ in range(int(t_end / dt)):
        u, v = update(u, v, D_u, D_v, mu_u, mu_v, k_1, k_2, k_3, dh, dt)
    return u

03-01. 拡散係数を変える¶

基本パラメータを固定し，活性化因子の拡散係数 D_u を 1.2〜1.42 の範囲で少しずつ変えてみる．

D_u_values = [1.2, 1.27, 1.3, 1.35, 1.37, 1.38, 1.4, 1.42]
fig, axes = plt.subplots(2, 4, figsize=(18, 9))
for ax, D_u in zip(axes.flat, D_u_values):
    ax.set_aspect("equal")
    ax.pcolormesh(x, y, simulate(D_u), cmap="gray")
    ax.set_title(f"D_u = {D_u:.2f}")
    ax.axis("off")
plt.show()

D_u を大きくしていくと，パターンの形が斑点（spots）から迷路状，そして縞（stripes）へ，さらに明暗が反転した穴状のパターンへと変化する．

03-02. パターンが生じる条件¶

チューリングパターン (Turing pattern)が生じるには，抑制因子が活性化因子より十分速く拡散する必要がある（ $D_v \gg D_u$ ）．

cases = [
    ("D_u=1.2, D_v=10 (pattern)", dict(D_u=1.2, D_v=10.0)),
    ("D_u=1.5 (activator too fast)", dict(D_u=1.5, D_v=10.0)),
    ("D_v=5 (inhibitor too slow)", dict(D_u=1.2, D_v=5.0)),
]
fig, axes = plt.subplots(1, 3, figsize=(14, 5))
for ax, (title, kw) in zip(axes, cases):
    ax.set_aspect("equal")
    ax.pcolormesh(x, y, simulate(**kw), cmap="gray")
    ax.set_title(title)
    ax.axis("off")
plt.show()

D_u を大きくしすぎる（活性化因子が速すぎる）か，D_v を小さくしすぎる（抑制因子が遅すぎる）と，パターンが生じない．

なお，初期のノイズはシードを固定していないので，実行するたびにパターンの配置は変わる．しかし斑点の大きさや間隔（パターンの特徴的なスケール）はあまり変わらない．これはパターンが初期条件ではなく，モデルのメカニズム（ここでは拡散係数の比）で決まるためである．

03-03. 生物の模様との対応¶

動物の体表模様（ヒョウの斑点，シマウマの縞など）の形成には，反応拡散系のようなメカニズムが関与していると考えられている．パラメータを変化させると斑点から縞へとパターンの型が変わり，これがヒョウ（斑点）とシマウマ（縞）のような模様の違いに対応しうる．

ギーラー-マインハルト系だけでも，パラメータによって斑点・縞・穴など多様なパターンが生じる．Gray-Scott モデルなど他の反応拡散系も，同様に多彩なパターンを示す．この章の課題では，パラメータや初期値を変えて様々なパタンを調べたり，特定の生物の模様の再現を試みたりする．

Solution to Exercise 1

例として D_u を変えると，斑点と縞という異なる型のパターンが得られる．

fig, axes = plt.subplots(1, 2, figsize=(10, 5))
for ax, D_u in zip(axes, [1.0, 1.3]):
    ax.set_aspect("equal")
    ax.pcolormesh(x, y, simulate(D_u), cmap="gray")
    ax.set_title(f"D_u = {D_u}")
    ax.axis("off")
plt.show()