反応拡散モデル - 数理生物学演習

拡散に反応項を加えると，ほぼ一様な状態から自発的に空間パターンが生じる．

import matplotlib.pyplot as plt
import numpy as np

02-01. ギーラー-マインハルト系¶

反応拡散系 (reaction-diffusion system)は，拡散する複数の因子が反応しながら空間パターンを形成する系である．ギーラー-マインハルト系では，活性化因子（activator） $u$ と抑制因子（inhibitor） $v$ の相互作用が次の式で表される．

\frac{\partial u}{\partial t} = D_u \nabla^2 u - \mu_u u + \frac{k_1 u^2}{v} + k_2

(1)

\frac{\partial v}{\partial t} = D_v \nabla^2 v - \mu_v v + k_3 u^2

(2)

右辺第1項は拡散，第2項は分解，残りが活性化や抑制を表す． $u$ は自己を活性化しつつ $v$ を増やし， $v$ は $u$ を抑制する．拡散の速い抑制因子（ $D_v \gg D_u$ ）が活性化を離れた場所で抑えることで，チューリングパターン (Turing pattern)が生じる．

各記号の意味は次のとおり．

記号	意味
`D_u`, `D_v`	拡散係数（活性化因子・抑制因子）
`mu_u`, `mu_v`	分解率
`k_1`	活性化因子の自己活性化（抑制因子で抑制）
`k_2`	活性化因子の基本生成
`k_3`	活性化因子による抑制因子の生産

02-02. 更新関数の実装¶

拡散項は前節と同じ周期境界条件のラプラシアンを使う． laplacian を再掲する．

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

反応項を加え， $u$ と $v$ を同時に更新する関数を定義する．モデルの式をそのままベクトル化して実装できる．

def update(u, v, D_u, D_v, mu_u, mu_v, k_1, k_2, k_3, dh, dt):
    """反応拡散モデル（ギーラー-マインハルト系）の1ステップ更新．

    周期境界条件を仮定している．

    Args:
        u, v: 活性化因子・抑制因子の場の配列
        D_u, D_v: 拡散係数
        mu_u, mu_v: 分解率
        k_1: 活性化因子の自己活性化（抑制因子で抑制）
        k_2: 活性化因子の基本生成
        k_3: 活性化因子による抑制因子の生産
    Returns:
        u_new, v_new: 次のステップの場の配列
    """
    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

02-03. シミュレーションの実行¶

パラメータを設定し，（1に微小なノイズを与えた）ほぼ一様な初期状態から時間発展させる．擬似乱数生成器は np.random.default_rng で作る．

D_u = 1.2  # 活性化因子の拡散係数（1〜1.5 くらいで動かす）
D_v = 10.0  # 抑制因子の拡散係数
mu_u = 1.0  # 活性化因子の分解率
mu_v = 1.0  # 抑制因子の分解率
k_1 = 1.0
k_2 = 0.05
k_3 = 1.0

dt = 0.02  # 時間の刻み幅
dh = 1.0  # 空間の刻み幅

x = np.arange(-25, 25, dh)
y = np.arange(-25, 25, dh)
xmesh, ymesh = np.meshgrid(x, y)

t_end = 1000.0  # t の終端（パターンが十分に発達するまで）
n_steps = int(t_end / dt)

# u, v を初期化
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)))

t_list = [0.0]  # 時刻を記録するリスト
u_list = [u.copy()]  # u を記録するリスト
v_list = [v.copy()]  # v を記録するリスト
for n in range(1, n_steps + 1):
    u, v = update(u, v, D_u, D_v, mu_u, mu_v, k_1, k_2, k_3, dh, dt)
    if n % 100 == 0:  # 100ステップごとに記録
        t_list.append(n * dt)
        u_list.append(u.copy())
        v_list.append(v.copy())

02-04. 可視化¶

最終状態の活性化因子 $u$ と抑制因子 $v$ を並べてプロットする．

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
for ax, field, name in zip(axes, [u_list[-1], v_list[-1]], ["u (activator)", "v (inhibitor)"]):
    ax.set_aspect("equal")
    pcm = ax.pcolormesh(xmesh, ymesh, field)
    fig.colorbar(pcm, ax=ax, shrink=0.75)
    ax.set_title(name)
fig.suptitle(f"t = {t_list[-1]:.0f}")
plt.show()

ほぼ一様だった初期状態から，斑点状のパターンが浮かび上がる．時間発展を並べると，パターンが形成されていく過程がわかる．

idx = np.linspace(0, len(u_list) - 1, 10).astype(int)
fig, axes = plt.subplots(2, 5, figsize=(16, 7))
for ax, k in zip(axes.flat, idx):
    ax.set_aspect("equal")
    ax.pcolormesh(xmesh, ymesh, u_list[k], cmap="gray")
    ax.set_title(f"t = {t_list[k]:.0f}")
    ax.axis("off")
fig.suptitle("Activator u over time")
plt.show()

Solution to Exercise 1

D_u を大きくすると斑点がつながり，縞状のパターンが出現する．

D_u = 1.39  # 活性化因子の拡散係数を大きくする

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 n in range(n_steps):
    u, v = update(u, v, D_u, D_v, mu_u, mu_v, k_1, k_2, k_3, dh, dt)

fig, ax = plt.subplots(dpi=120)
ax.set_aspect("equal")
pcm = ax.pcolormesh(xmesh, ymesh, u)
fig.colorbar(pcm, ax=ax, shrink=0.8)
ax.set_title(f"D_u = {D_u}")
plt.show()