現象の数理モデル構築法
長山 雅晴(北海道大学 電子科学研究所 附属社会創造数学研究センター)∗

1. はじめに現象の数理モデリングとは何でしょうか？数理モデリングとは，実験や観察から示唆される法則や仮定から，対象とする現象を数学の問題として定式化することではないでしょうか．数理モデルが現象を定性的に再現できれば，数理モデルを解析することで，現象の起因となる原理を理論から示唆することが可能となります．それだけでは，現象を理解したことにならないことは十分承知していますが，支配原理が不明な現象に対して支配原理を見つける一端になることが期待できると考えられます．この原稿では時間変化していく現象の数学モデルの構成法について説明します．残念ながら，支配原理が不明瞭な現象に対する確立した数理モデリング手法はないので，私が数理モデルを構築するために考えている方法を説明したいと思います．時間変化していく現象を記述する数学の道具として
セル・オートマトン，漸化式，微分方程式
等があります．では，これらの道具を使ってどのようにして数理モデルを構築すればよいのでしょうか？数理モデルを構成するための絶対的な原理というものはないと思いますが，これまでの経験から次の基本原理によって数理モデルは構成されるように私は考えています：

数理モデリングの基本原理
単位時間あたりのある（ものの）量の変化量 = 流入量 − 流出量= 生成量 − 消費量

ここで，ある（ものの）量とは化学物質濃度，物質密度，個体群密度，温度，速度などを指します．時間とともに量が変化（状態が遷移）する現象の数理モデルはこの基本原理によって記述されていると思ってよいでしょう．ここでは微分方程式を用いた
数理モデリングについて解説していきたいと思います．このとき数理モデリングの基
本原理の左辺はつぎのように記述できます：


ある量 x の単位時間あたりの変化量 =


dx
dt





あとは対象とする現象の流入量と流出量をあたえることができれば，数理モデルの原型をつくることができます．重要なポイントは「流入する（もの）の量」と「流出する（もの）の量」を見つけることになります．そのためには，対象とする現象に対して自然な形で仮説を立て，その仮説から右辺を与えなければなりません．
簡単な例としてある 1 種類の生物の個体数の増減を表すモデルを考えてみましょう．
このモデルは次のように考えることができます：∗〒 060-0807 札幌市北区北 12 条西 7 丁目 中央キャンパス総合研究棟 2 号館
e-mail: nagayama@es.hokudai-u.ac.jp

単一個体群の増減モデル
単位時間あたりの個体群の変化量 = 出生数 − 死亡数

ここで，時刻 t における個体群密度を x(t)，単位時間あたりの出生率と死亡率をそれぞれ a, b（a, b は定数）とすると，上記の数理モデルは次のような微分方程式で書くことができます：
dx= ax − bx.
(1)dt
このモデル方程式は，マルサスの人口増減モデルと呼ばれています [1]．モデル方程式(1) において，ある時刻 t0 での個体群密度を x0 とすると，(1) の解はx(t) = x0 exp((a − b)(t − t0 ))となります．従って


 +∞, a − b > 0,
lim x(t) =x0 ,a − b = 0,
t→+∞

 0,
a−b<0

(2)

となることがわかります．すなわち，出生率より死亡率が低い場合，個体群密度は増加し，出生率より死亡率が高い場合は個体群密度は減少することがわかります．この
例のように大学 1 年生で習う簡単な線形常微分方程式も数理モデル 1 となり，そこから個体群密度の時間変化が予測されるのです．ここで主張したいことの一つは，与えられた微分方程式の変数や定数が現象の言葉として意味を持てば，その微分方程式は数理モデルになるのだということです．
数理モデル (1) において，出生率より死亡率が低い場合，(2) から個体群密度は抑制されることなく増加し，発散していくため，個体群増減モデルとしては不十分である
ことがわかります．そのため，個体群密度が高くなると出生率が低下するという仮定を加えた数理モデル
(
x)
dx
=a 1−
x.
(3)
dt
c
が提唱されています．このモデルはロジステック方程式と呼ばれています [1]．ここで，c は環境収容量と呼ばれる定数です．
この解説書では，主に化学反応に対する数理モデリングについて説明していきます．
その応用として燃焼合成反応の数理モデリングについて解説します．最後に，化学反応の数理モデル化原理と同じ「質量作用の法則」を用いて書き下すことかできる生物種に対する個体群密度の数理モデリングについて説明したいと思います．

2. 化学反応の数理モデル
この節ではまず化学反応に対する数理モデリングについて説明し，次にその応用として，Belousov-Zhabotinsky (BZ) 反応と個体群動態の数理モデリングについて述べたいと思います．化学反応の数理モデルは反応の素過程から構成することができます．素過程とは化学反応が完了するまでに起こる一つ一つの反応過程を記述している化学反応のことです．例えば，次のような水素分子と酸素分子から水が生成される化学反応式k

2H2 + O2 −→ 2H2 O
1

ここで構成した数理モデルが適切なモデルであるかないかは別ですが…

に対する簡略化された素過程 2 の反応式は{
k1
H+ + O2− −→ OH−
2
H+ + OH− −
→
H2 O

k

(4)

となり，化学反応式とその素過程の反応式は異なります．このような化学反応式と素過程の反応式が異なるような反応は複合反応と呼ばれています．一方，化学反応式と素過程の反応式が同じような反応は単純反応と呼ばれています．では，素過程の反応
式からどのように数理モデルを構築すればよいのでしょうか．次のような素過程の反応式を仮定しましょう：kma A + mb B −→ mp P + mq Q
(5)
ここで，ma , mb , mp , mq は化学量論係数 3 ，k は反応速度定数と呼ばれています．化学
反応では，化学量論比を一定に保ちつつ「質量作用の法則」に従って反応が進むことが知られているため，その原理に基づいて数理モデリングを行います．反応式 (5) において化学種 A, B は反応によって消費されるので，それらをまとめると数理モデルは次のように記述されます：
1 da
1 db
=
= −kana bnb .
ma dt
mb dt
ただし，a, b はそれぞれ化学種 A, B の濃度を，na , nb は反応次数を表します．反応式 (5)において化学種 P, Q は A と B の反応によって生成されるので，それらをまとめると数理モデルは次のように記述されます：

1 dp
1 dq
=
= kana bnb .
mp dt
mq dt
ここで，p, q は化学種 P, Q の濃度を表します．一般に反応のべき係数 na , nb は化学量論
係数 ma , mb と一致しない場合もありますが，以下の記述では化学量論係数と反応のべき係数は同じであると仮定します 4 ．ここまでの説明で最初に説明した素過程の反応式(4) の数理モデルを構成することができるので挑戦してみてください．
次に，可逆反応に対する数理モデルを考えます．ここでは，次のような素過程の反応式を考えます：k+ma A + mb B ⇆ mp P

(6)

k−

この反応では，化学種 A は B と反応して消費され，化学種 P の自己分解から生成されるので，反応式 (6) に対する化学種 A の数理モデルはda
1 da= −k+ ama bmb + k− pmp から= ma (−k+ ama bmb + k− pmp )ma dtdtと記述できます．化学種 B も同様に記述できます．一方，化学種 P は A と B の反応により生成され，P 自身の自己分解反応で減少するので，反応式 (6) に対する化学種 P の2

実際の素過程はイオン化のプロセスや同時に進行する複数の反応経路を含むため複雑になります．
反応に必要な mol 数と考えてもいいです．すなわち，化学種 A が ma mol と化学種 B が mb mol から
化学種 P が mp mol と化学種 Q が mq mol できると考えられます．
4
素過程が正しく記述されていても，数理モデルが化学反応の定性的性質を再現できないようであれば，化学量論係数と反応のべき係数は異なっている可能性があります3

数理モデルは

1 dp
dp= k+ ama bmb − k− pmp から= mp (k+ ama bmb − k− pmp )mp dtdtと記述することができます．このように，化学反応に対する素過程の反応式が与えられれば，数理モデルは比較的簡単に記述できることがわかります．ここで示した化学反応の例では，数理モデルの濃度 a(t0 ), b(t0 ), p(t0 ), q(t0 ) を与えると，初期濃度に応じた平衡濃度に収束していくことが数理モデルの解析から確かめることができます．
最後に，自己触媒反応の数理モデルについて説明します．自己触媒反応とは，例えばkma A + B −→ (ma + m)A(7)のように化学種 A 自体が触媒となって A を増加させる反応です．この反応の数理モデ
リングは次のように考えます．まず，(7) の左辺の反応によって化学種 A は減少することから，左辺の反応から得れる化学種 A の濃度変化を表す数理モデルは

1 dada= −kama b から
= −ma kama b
ma dt
dtと記述できます．一方，(7) の左辺の反応によって化学種 A は増加することから，右辺の生成過程から得られる化学種 A の濃度変化を表す数理モデルはdada1= kama b から= (ma + m)kama b
ma + m dt
dtとなります．従って，全体として化学種 A の時間変化を表す数理モデルは「数理モデリングの基本原理」を用いると

da= −ma kama b + (ma + m)kama bdtとなるので，最終的に

da= mkama b
dtと記述できます．つまり，化学反応式の右辺の化学量論係数から左辺の化学量論係数を引いた差だけ増加率が変わったモデル方程式が導出されるのです 5 ．反応式の中に自己触媒反応が含まれている場合は，単純に平衡濃度に収束しないことがあり得ます．次節はこの自己触媒反応が含まれていることに起因して振動する化学反応の数理モデルについて紹介します．
2.1. ベローソフ・ジャボチンスキー反応の数理モデリング
この節では，前節で説明した素過程の反応式を用いた数理モデリングを応用して，BZ
反応 [2, 3, 4] の数理モデルを構成してみましょう．BZ 反応は酸化と還元を交互に起こし続ける振動反応の例して有名な酸化還元反応です．BZ 反応の機構は FKN メカニズ
ム [5] とよばれており，10 個の素過程による振動メカ二ズムが提案されています．そこ5

実際の化学反応では，ほとんどの場合 m = 1 となります．

から，可逆反応を無視してさらに単純化した次の 5 つの素過程が提案されました [6]：

k1


A+Y −
→X+P



k

2
X + Y −
→
2P


k3
A+X−
→ 2X + 2Z (自己触媒反応)


k
2X −
4

→
A+P




k
5
Z + B −→hY

BZ 反応では化学種 A と B は反応によって濃度変化がほとんど見られないことから A と
B の濃度は定数と仮定します．このとき，化学種 X, Y, Z の濃度を x(t), y(t), z(t) として，化学種 A,B の濃度を a, b(a, b は定数) とすると，BZ 反応のモデル方程式はdx= k1 ay − k2 xy + k3 ax − 2k4 x2 ,


dt



dy
(8)
= −k1 ay − k2 xy + k5 zbh,

dt





 dz = 2k3 ax − k5 bz
dtと記述されます. この数理モデルは，反応過程を提案した 3 名の所属大学の名前をとってオレゴネーターモデルと呼ばれています．数理モデル (8) から縮約（特異摂動法）することによって，よく知られている 2 変数の BZ 反応モデル

dx
q−x


= fz
+ x(1 − x),
ε
dt
x+q
(9)

dz

 =x−zdtが導出されています [2, 7]．2 変数モデル (9) は Keener-Tyson モデルと呼ばれており，(8) に対して無次元化（単位の規格化）を行い，特異摂動法によって変数 y を断熱的に消去することによって導出することができます [2]．数理モデルを数学として取り扱う
ためには無次元化は必須であり，実験から示唆されるデータからモデル方程式を縮約してより簡単な数理モデルを導出することは重要な作業なので，以下では BZ 反応モデ
ルを例に取り，それらの手続きを説明したいと思います．

2.1.1. 無次元化
無次元化とはモデル方程式の独立変数および従属変数が持つ単位をすべて 1 に統一するような変数変換のことです．(8) が数理モデルとして成立していると仮定すると両辺は同じ単位を持たないといけません．左辺は単位時間 [s] 当たりの濃度変化 [mol/l3 ] となっていますので、右辺も同じ単位を持たないといけません．そのためには，ki の単
位は [1/(s·(mol/l3 )] となることがわかります．そこで次のような変数変換を行います：

τ = k5 bt, x =k4k2x, y =
y,
k3 a
k3 a

z=

k4 k5 bz.
2(k3 a)2このとき，数理モデル (8) は x → x, y → y, z → z と書き換えるとdx
= qy − xy − +x(1 − x),
ε


dτ



dy
ε′
= (−qy − xy + f z),

dτ





 dz = x − zdτとなります．ただし、
k4 k5 bk1 k4
k5 b
, ε′ =, q=, f = 2hε=
k3 a
k2 k3 a
k2 k3

(10)

です．(10) は変数およびパラメーターの単位はすべて 1 となっています．これで無次元化が完了しました．

2.1.2. 特異摂動法
ここでは (10) に対して特異摂動法によって変数 y を断熱的に消去することによって (10)k4 k3を導出する方法を説明します [2]．q = 2k2 ( )( ) = 8 × 10−4 となることが，実験か
k1 k4
らわかっています．実験では反応速度の測定は困難であっても反応速度比は測定可能な場合があり，BZ 反応でも反応速度比を計測することができて，0 < ε′ ≪ ε < 1 を仮定することができます．(10) から ε′ → 0 とすると第 2 式から

fzx+qとなり，これを 1 式と 3 式に代入すると，(9) が導出されます．このような手続きを物理では y を断熱消去をするということがあります．ここでは，第 2 式で ε′ → 0 とすると微分方程式の最高階微分の係数を 0 としていることから特異摂動と呼ばれる作業を行っています．このような手続きを形式的に行ってもよいことを説明したいと思います．モデル方程式 (10) に対する初期値を (x, y, z)(0) = (x0 , y0 , z0 ) とします．最初に，非常に短い時間 [0, ε′ ] の間に (10) の解の振る舞いを調べます．そのためにtτ′ = ′
ε′と変数変換します．これによって，[0, ε ] の短い時間での解の変化を調べることができます．このとき (10) はdxε ′ = ε′ (qy − xy + x(1 − x)),


dτ



dy
(11)
= −qy − xy + f z,
′

dτ





 dz = ε′ (x − z)dτ ′となります．ここで，ε′ → 0 とすると
dxε ′ = 0,


dτ



dy
(12)
= −qy − xy + f z,
′

dτ





 dz = 0dτ ′y=を得ます．従って，第 1 式と第 3 式からx(τ ) ≡ x0 ,z(τ ) ≡ z0となり，この定数解を第 2 式に代入すると

dy= −qy − x0 y + f z0
dτ ′

(13)

t
から limτ ′ = +∞ となるので，(13) の漸近挙動を調べればよいこと
′
′
ε →0
ε
がわかります．(13) の平衡解は
f z0
y=x0 + qを得ます．τ ′ =であり，この平衡解は漸近安定であることが容易にわかるので、任意の初期値 y(0) 対
して
f z0
′
lim
y(τ
)
=
τ ′ →+∞
x0 + qとなることがわかります．従って，y(τ ) は短い時間 [0, ε′ ] でy(τ ) =f z(τ )x(τ ) + q

(14)に近づくことがわかります．このことは，短い時間 [0, ε′ ] 以降の時間では y(τ ) が x(τ ), z(τ )によって決まることを示唆しています．従って，(10) において ε′ → 0 としたのち，(14)を (10) に代入することで 2 変数 Belousov-Zhabotinsky 反応モデル

dx
q−x


= fz
+ x(1 − x) ≡ F (u, z),
ε
dτ
x+q
(15)

dz


= x − z ≡ G(u, z).
dτを得ることができます．このモデルの右辺のヌルクラインを図に表しておきます．パ
ラメーター f, q 値に依存して平衡点の安定性が変わることがわかります．

2.2. 燃焼合成反応の数理モデリング
この節では，燃焼合成とよばれる材料合成の簡単な数理モデリングについて解説します．例としてチタンとカーボンの混合粉からチタンカーボンを生成する反応系

Ti + C → TiCを考えます．混合粉を円柱状にかため（圧粉体とよばれる），その上端をカーボンヒー
トで熱すると，約 1800 度で化学反応が始まります．この反応は発熱反応系であり，反応熱によっては 2200 度以上に上昇し，圧粉体の上端から下端まで反応熱の拡散により未反応部分が順次反応を起こします．そして自動的に反応が終了し、TiC の化合物ができます．この反応をプロセスを記述する数理モデルをつくってみましょう！ここではk(T )mA + nB −−→ P + 熱

0.14

v

G=0

u
0.1

0.06

F=0
0.02
0
0

0.2

0. 4

0. 6

0.8

1

図 1: (15) のヌルクライン, f = 2.0, q = 0.02 のとき，平衡点は安定平衡点となっている.
0.14

v
u

0.1

F=0

0.06

G=0

0.02
0
0

0.2

0.4

0. 6

0.8

1

図 2: (15) のヌルクライン, f = 3.0, q = 0.02 のとき，平衡点は安定平衡点となっている.
という一般的な発熱反応を考えます．材料合成からの要請として「完全に反応が完了した時に A と B はすべて P になっている」ことを条件とします．a を物質 [A] の密度，bを物質 [B] の密度，T :温度としたときモデル方程式は


da


= −mam bn k(T ),


dt



db
(16)
= −nam bn K(T ),

dt





 dT = αam bn K(T )dtと記述できます．ここで，α は熱変換率，K(T ) はアレニウスの反応率とよばれておりk exp(− Eapp ), T ≥ T ,0ig
RT
K(T ) =
(17)
0,
T <T
ig

と記述されます．ここで，Tig は反応開始温度，Eapp は見かけの活生化エネルギー，k0は指数前定数です．このとき (16) にたいして，n ×第一式‐m × 第二式とすると

n

da
db
−m =0
dt
dt

となる．ここで，両辺を t で積分すると
∫ t
da
db
(n − m )dt = 0
dt
dt
0

となり

na(t) − mb(t) = na(0) − mb(0)(18)

を得ます．反応終了時にはすべて生成物になる条件から

lim (na(t) − mb(t)) = 0t→∞が必要です．そのためには

na(0) = mb(0)が成立していればよいので，初期条件として

a(0) : b(0) = m : nを満たす必要があります．すなわち初期条件においてモル比を a(0) : b(0) = m : n で圧粉体を制作すればよいことがわかります．合成反応の初期条件として

a(0) : b(0) = m : nとおくと，(18) からna(t) = mb(t)を得ます．従ってna(t)mとなり，(16) の第 1 式と第 2 式は同じ方程式
b(t) =dan= −mam ( )n an k(T )dtmとなるため，モデル方程式は

da
n


= −m( )n am+n k(T ),

dt
m

n
dT


= α( )n am+n k(T )
dt
m

(19)

と 2 変数の数理モデルに縮約されます．ここで圧粉体での熱伝導を考えると，モデル方
程式は

da
n


= −m( )n am+n k(T ),

dt
m
(20)

dT
n n m+n= dT ∆T + α( ) a
k(T )
dt
m
と記述されます．ここで dT は圧粉体での熱伝導係数です．圧粉体中での粒子の拡散速度スケールと反応速度のスケールを考えると，粒子の拡散は無視ほど小さいため，反応物の拡散は導入しないことにします．圧粉体は円柱状のため円柱領域で計算することを考えます．円柱領域 Ω での数値計算を行うために，直交座標系 (x, y, z) を円筒座標系 (r, θ, z) に変換します．この変数変換は

x = r cos θ, y = r sin θ, z = z(21)

として与えられます．このとき，円柱領域 Ω は Ω = {(r, θ, z|0 < r < R, 0 < θ < 2π, 0 <z < zL } と表現でき

∆=

∂2
∂2
∂2
∂2
1 ∂
∂2
∂2
+
+
=
+
+
+
(= ∆r,θ,z )∂x2 ∂y 2 ∂z 2
∂r2 r ∂r ∂θ2 ∂z 2となります．このときモデル方程式 (20) は

da
n


= −m( )n am+n k(T ),

dt
m

dT
n


= dT ∆r,θ,z T + α( )n am+n k(T )
dt
m

(22)

と表現できます．円柱の上端と下端について熱入出がないと仮定して，斉次ノイマン
境界条件
∂T
∂T
(r, θ, 0) = 0,(r, θ, L) = 0, 0 < r < R, 0 < θ < 2π,∂z∂zを与えます．円柱側面での境界条件は熱のロスが燃焼合成に与える影響を調べるために，物理的に自然な条件として• アルゴン (不活性気体) 気体中での熱対流による放熱∂T(R, θ, z) = λ(Tθ − T ), 0 < θ < 2π, 0 < z < zL ,
∂r

(23)

• 真空中における輻射による放熱∂T(R, θ, z) = εb(T 4 − T04 ), 0 < θ < 2π, 0 < z < zL
∂r

(24)

の 2 つの境界条件を考えます．また，燃焼合成反応における物理的に制御可能なパラメーターは Eapp（見かけの活生化エネルギー），R（円柱の半径），T0（外気温）です．
これらのパラメーターを変化させることで燃焼合成反応の反応面がどのような形態を取るか調べていくことが，重要な課題となります．円筒座標系での数値計算法については私の Web ページにある反応拡散系数値計算入門をご覧下さい [13]．

3. 個体群動態の数理モデリング
数理生態学で用いられている個体群密度の数理モデルの構成では質量作用の法則を用いており，化学反応モデルの構成法と類似しているので，化学反応モデルの応用として取り上げたいと思います．ここで導出される数理モデルの解説書は多く出ているので，そちらも是非参考にしてみてください（例えば，[8, 9] です）．

3.1. ロジスティックモデル
ある個体種 X が，個体種 N を捕食し増加する問題を考えます．ここでは 1 個の X が 1 個の N を捕食して増加すると仮定し，単位時間あたりの捕食率を α，N を 1 個体を捕食したときに X が増加する割合を β とします．このとき化学反応式を応用すると
αX+N−→ (1 + β)X

と書くことができます．質量作用の法則を適用して数理モデリングを行うとdx

= αβxn,dt dn = −αxn
dt

(25)

と書くことができます．ここで x, n は個体種 X，N の個体群密度を表しています．この数理モデルは最も簡単な捕食者・非捕食者モデルとなっています．このとき，X 種，N種の初期個体群密度 x(0), n(0) をそれぞれ x0 , n0 とすると，(25) の 2 式から

d(x + βn) = 0dtを得ることができx(t) + βn(t) = x0 + βn0となることがわかります．この式から n(t) を (25) の第一式に代入すると()
1dx= α(x0 + βn0 ) 1 −x x
dt
x0 + βn0(26)となり，1 節で示したロジステック方程式 (3) が導出されます．(26) から x0 > 0, n0 > 0
であれば
lim x(t) = x0 + βn0t→∞となることがわかり，limt→∞ n(t) = 0 も導かれます．
数学の問題では x0 や n0 の符号に合わせて場合分けをすればよいのですが，個体群密度のモデル方程式なので x(t) や n(t) は 0 以上の値を取らなければ密度としての意味をもたなくなります．従って，初期値も 0 以上でなければ意味がありません．また，初期値を正の値に取ったとしても，ある時刻で個体群密度の値が負になるような数理モデ
ルであれば，その数理モデルは破綻しているといえるので，数理モデルを構成する際
は解の正値性も調べる必要があります．ところで，数理モデル (25) と (26) に解析結果
から，有限の餌を捕食しながら増殖する捕食種 X はその個体群密度に上限があることがわかるので，モデル方程式 (3) を作るときに導入した「個体群密度が高くなると出生率が低下する」という仮定は，餌が有限であることと同値であることがわかります．

3.2. 2 種競争系モデル
ある 2 種の個体種 X，Y が，同じ個体種 N を捕食し増加する問題を考えてみましょう．
ここでは 1 個の X や Y が 1 個の N を捕食して増加すると仮定し，X，Y の単位時間あたりの捕食率を αx , αy とし，N を 1 個体を捕食したときに X，Y が増加する割合を βx , βy
とすると

αx

X + N −→(1 + βx )XαyY + N −→ (1 + βy )Yという式が得られます． このとき数理モデルはdx= αx βx xn,
dt

(27)

dy
= αy βy yn,
dt

(28)

dn
= −αx nx − αy nx
dt

(29)

となります．ここで x, y, n は個体種 X，Y，N の密度を表しています．このとき，X 種，
Y 種，N 種の初期個体群密度 x(0), y(0), n(0) をそれぞれ x0 , y0 , n0 とします．(27)–(29)からd(βy x + βx y + βx βy n) = 0
dtとなり(βy x + βx y + βx βy n)(t) = βy x0 + βx y0 + βx βy n0(30)を得ます．(30) を (27) と (28) に代入して整理すると

dx

= ax (1 − bx x)x − kx xy,dt dy = a (1 − b y)y − k xy
y
y
y
dt

(31)

となります．ただし

ax = (βy x0 + βx y0 + βx βy n0 )

αx
βy
βx
, bx =, kx = αx ,βyβy x0 + βx y0 + βx βy n0
βyay = (βy x0 + βx y0 + βx βy n0 )

αy
βx
βy
, by =, ky = αy .
βx
βy x0 + βx y0 + βx βy n0
βx

モデル方程式 (31) は 2 種競争系モデルと同じ方程式になります．2 種競争系モデルに現れる種間競争項 xy は種 X と種 Y が直接競合する仮定から与えられていますが，種 Xと種 Y が同じ種 N を捕食するという仮定からも種間競争項を導出することができます．
この結果から，餌を取り合うという競争は，餌である被捕食種 N の動態を微分方程式として表すことなく，種間競争として表現できることがわかります 6 ．ここで説明した
数理モデリングを 1 種の被捕食種 N を k 種の捕食種が捕食し増殖する問題に拡張すると，最後には k 種の競争系モデルが導出できるので，是非挑戦してみてください．

4. まとめ
本稿では，数理モデリングを行うための基本原理を説明し，その原理に基づいて非常にシンプルな仮定での数理モデリングの例を説明し，化学反応モデル，個体群動態モデルの構成について説明しました．本稿では説明しませんでしたが，個体群動態のほかに感染症モデルも化学反応モデルと似た形式で記述できます [9]．また，生体内で生じる生化学反応（例えば，酵素分解反応）や細胞内でのカルシウムイオン興奮なども6

ただしこの方法で構成された 2 種競争系モデルの競争係数は対称な形になるため，2 種間での競争の強弱をつけることができません．もっと改良した導出方法を考える必要があります．

同じような方法で導出することができます [10]．ここで説明した数理モデル以外に，ろうそく火炎振動子の同期現象 [11] やメトロノームの同期現象に対する数理モデリング
も物理法則と化学反応の数理モデルを応用すると構築できます．そして，最近研究している表皮の数理モデリング [12] も同様に構築することが出来ますが，ここでは割愛
させて頂きます。これらの数理モデリング構築について詳しく知りたい方は私の Webページから PDF ファイルをダウンロードしてください [14]．最後になりますが，2 つのことを強調しておきたいと思います．1 つは，与えられた微分方程式の変数や定数が現象の言葉として意味を持てば，その微分方程式は数理モデルになるということです．
微分方程式の研究が何の役に立つの？とよく聞かれますが，研究対象とする微分方程式が現象としての意味を持てばその答えは明白だと思います．もう 1 つは，数理モデリングとは数式による現象の模倣であるということです．基本原理がわからない現象に対して，基本原理の一端を理論的に捉えるための道具なのです．基本的に仮説を与える時点で主観が入るため，どうしても数理モデリングは客観性を失います．従って，数理モデリングによって現象を再現できたからといって現象の数理モデルができたとは言えないのです．与えた仮説から可能な限り主観を取り除く努力をすることが必要
なので，数理モデルから示唆されたメカニズムや予測を実験で検証することが重要になってきます．そのためには数値計算をする必要があります．数値計算法についての簡単な解説はにありますので参考にして下さい [13]．

参考文献[1] “微分方程式で数理モデルをつくろう”，デヴィッド・バージェス，モラル・ボリー，日本評
論社[2] “非平衡系の科学 III”，三池 秀敏，森 義仁，山口 智彦，講談社.
[3] B. P. Belousov, in Collections of abstracts on radiation medicine(Medgiz, Moscow), 1958,
145–.
[4] A. M. Zhabotinsky, Bioﬁzika, 9, 1964, 306–.
[5] R. J. Field, E. Körös and R. M. Noyes, Oscillations in chemical systems. II. Thorough analysis of temporal oscillation in the bromate-cerium-malonic acid system,
J. Am. Chem., 94(25), 1972, 8649–8664.
[6] R. J. Field and R. M. Noyes, Oscillations in chemical systems. IV. Limit cycle behavior
in a model of a real chemical reaction, J. Chem. phys., 60(5), 1974, 1877–1884[7] J. P. Keener and J. J. Tyson, Spiral wave in the Belousov-Zhabotinskii reaction, Physica
D 21, 1986, 307–324.
[8] “数理生物学（個体群動態の数理モデリング入門）”, 瀬野裕美，共立出版[9] “マレー数理生物学入門”，James D. Murray，丸善出版[10] “数理生理学 上，下”，J. Keener，J. Sneyd，日本評論社[11] H. Kitahata, J. Taguchi, M. Nagayama, T. Sakurai, Y. Ikura, A. Osa, Y. Sumino,
M. Tanaka, E. Yokoyama and H. Miike, Oscillation and synchronization in combustionof candles: Radiation coupling could induce synchronization of oscillatory combustion,
Journal of Physical Chemistry (A), 113(29), 2009, 8164–8168.
[12] Y. Kobayashi, Y. Sawabu, H. Kitahata, M. Denda, M. Nagayama, Mathematical model
for calcium-assisted epidermal homeostasis, Journal of Theoretical Biology, 397, 2016,
52–60.
[13] http://mmc01.es.hokudai.ac.jp/wp/wp-content/uploads/2019/11/Intro Numeric RD.pdf[14] http://mmc01.es.hokudai.ac.jp/wp/wp-content/uploads/2019/11/Math Modeling Tutorial.pdf