書記が数学やるだけ#902 一成分の反応拡散方程式の例

前回扱った神経細胞のモデルの延長として，一成分の反応拡散方程式の例を見ていく。

問題

説明

拡散方程式とは、物質やエネルギー、情報などが空間的に広がっていく過程を記述する偏微分方程式であり、熱伝導・粒子拡散・化学反応・生物個体群の空間分布など、さまざまな自然現象の基本モデルとして広く用いられている。

数理モデルにおける反応拡散方程式の例として、フィッシャー＝コルモゴロフ方程式やアレン＝カーン＝南雲方程式がある。これらはいずれも、空間拡散と局所反応の競合によって、時空間的なパターン形成や進行波を記述する代表的な方程式である。フィッシャー＝コルモゴロフ方程式は、生物個体群の拡散と増殖を表すモデルとして導入されたものであり、アレン＝カーン＝南雲方程式は、相転移や神経興奮などの二値的な状態変化を記述するモデルとして知られている。

神経伝達の過程をモデル化することで、反応拡散方程式の一種であるケーブル方程式に帰着できる。ケーブル方程式は、神経細胞の軸索や樹状突起を電気的な伝導路（ケーブル）として捉え、膜電位の時間変化と空間拡散を記述するものである。

解答

ケーブル方程式は、電荷の変化および膜電位の変化に伴う電流の出入りを考慮し、それらの関係式を整理することで導出できる。

反応拡散方程式の例として、フィッシャー＝コルモゴロフ方程式を考える。この方程式は、式を無次元化し、進行波解を仮定したうえで線形化することで、線形微分方程式の形に変形できる。さらに、この線形化された方程式の解が存在する条件から、波速の最小値が決定される。

相平面解析により、平衡点として安定ノード (0,0) と サドル点 (1,0) が存在することが示される。

このコードは、フィッシャー＝コルモゴロフ方程式における進行波解の数値的形状を可視化するもので、波速 c の違いによる形状の変化に注目している。計算結果から、波速が大きくなるほど波の前線が急峻になり、前線の傾きや広がりが変化することが確認できる。特に最小波速 cmin =2 においては、最も緩やかな前線が形成されることが分かる。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_bvp

# 1. 常微分方程式: U'' + (c/D) U' + (r/D) U(1-U) = 0
#    ここでは無次元化して D = r = 1 とする
def f(z, U, c):
    # U[0] = U, U[1] = U'
    dUdz = np.zeros_like(U)
    dUdz[0] = U[1]
    dUdz[1] = -c * U[1] - U[0] * (1 - U[0])
    return dUdz

# 境界条件: U(-L) = 1, U(L) = 0
def bc(Ua, Ub):
    return np.array([Ua[0] - 1, Ub[0]])

# 計算領域
L = 20
z = np.linspace(-L, L, 1000)

# 波速リスト
c_list = [2, 2.5, 3, 3.5, 4]

plt.figure(figsize=(8,5))

for c in c_list:
    # 初期推定
    U_init = np.zeros((2, z.size))
    U_init[0] = 0.5*(1 - np.tanh(z/5))  # 初期形状の推定
    # solve_bvpでBVPを解く
    sol = solve_bvp(lambda z,U: f(z,U,c), bc, z, U_init)
    plt.plot(sol.x, sol.y[0], label=f'c = {c}')

plt.xlabel('z = x - ct')
plt.ylabel('U(z)')
plt.title('Fisher–KPP')
plt.grid(True)
plt.legend()
plt.show()

アレン＝カーン＝南雲方程式は、フィッツフュー＝南雲モデルの簡略化の一例である。この方程式に進行波解を仮定することで、ハクスリー解と呼ばれる単調な波形解の存在が示される。ハクスリー解は、二値的な状態変化や相転移現象を表すモデルにおいて、安定な前線の伝播を記述する代表的な進行波解である。

本記事のもくじはこちら：