現象の数理モデリング〜 化学反応モデルを中心に 〜

長山 雅晴(北海道大学 電子科学研究所 附属社会創造数学研究センター)

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

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

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

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



あとは対象とする現象の流入量と流出量をあたえることができれば，数理モデルの原型をつくることができます．重要なポイントは「流入する（もの）の量」と「流出する（もの）の量」を見つけることになります．そのためには，対象とする現象に対して自然な形で仮説を立て，その仮説から右辺を与えなければなりません．
簡単な例としてある 1 種類の生物の個体数の増減を表すモデルを考えてみましょう．
このモデルは次のように考えることができます：2018 年度秋季総合分科会市民講演会（9 月 23 日）






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

ここで，時刻 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) から個体群密度は抑制されることなく増加し，発散していくため，個体群増減モデルとしては不十分であることがわかります．そのため，個体群密度が高くなると出生率が低下するという仮定を加えた数理モデル
dxx= a(1 − )x.
(3)
dtcが提唱されています．このモデルはロジステック方程式と呼ばれています [1]．ここで，c は環境収容量と呼ばれる定数です．

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

k12(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 dp1 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 の数理モデルは

da1 da= −k+ ama bmb + k− pmp から= ma (−k+ ama bmb + k− pmp )ma dtdtと記述できます．化学種 B も同様に記述できます．一方，化学種 P は A と B の反応により生成され，P 自身の自己分解反応で減少するので，反応式 (6) に対する化学種 P の数理モデルはdp
1 dp= k+ ama bmb − k− pmp から= mp (k+ ama bmb − k− pmp )mp dtdtと記述することができます．このように，化学反応に対する素過程の反応式が与えられれば，数理モデルは比較的簡単に記述できることがわかります．ここで示した化学3反応に必要な mol 数と考えてもいいです．すなわち，化学種 A が ma mol と化学種 B が mb mol から
化学種 P が mp mol と化学種 Q が mq mol できると考えられます．
4
素過程が正しく記述されていても，数理モデルが化学反応の定性的性質を再現できないようであれば，化学量論係数と反応のべき係数は異なっている可能性があります反応の例では，数理モデルの濃度 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 bma dtdtと記述できます．一方，(7) の左辺の反応によって化学種 A は増加することから，右辺の生成過程から得られる化学種 A の濃度変化を表す数理モデルは1da
da= kama b から= (ma + m)kama bma + m dtdtとなります．従って，全体として化学種 A の時間変化を表す数理モデルは「数理モデリングの基本原理」を用いると

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

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

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

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 bzdtと記述されます. この数理モデルは，反応過程を提案した 3 名の所属大学の名前をとってオレゴネーターモデルと呼ばれています．数理モデル (8) から縮約（特異摂動法）することによって，よく知られている 2 変数の BZ 反応モデル
dxq−x= fz+ x(1 − x),εdtx+q(9)dz =x−zdtが導出されています [2, 7]．2 変数モデル (9) は Keener-Tyson モデルと呼ばれており，(8) に対して無次元化（規格化）を行い，特異摂動法によって変数 y を断熱的に消去することによって導出することができますが，ここではその詳細は省略します. 興味がある方は [2] を参照してください．

2.2. 個体群動態の数理モデリング
数理生態学で用いられている個体群密度の数理モデルの構成では質量作用の法則を用いており，化学反応モデルの構成法と類似しているので，化学反応モデルの応用として取り上げたいと思います．ここで導出される数理モデルの解説書は多く出ているので，そちらも是非参考にしてみてください（例えば，[8, 9] です）．ある個体種 X が，個体種 N を捕食し増加する問題を考えます．ここでは 1 個の X が 1 個の N を捕食して増加すると仮定し，単位時間あたりの捕食率を α，N を 1 個体を捕食したときに X が増加する割合を β とします．このとき化学反応式を応用すると
αX+N−→ (1 + β)Xと書くことができます．質量作用の法則を適用して数理モデリングを行うと
dx

= αβxn,dt dn = −αxndt(10)と書くことができます．ここで x, n は個体種 X，N の個体群密度を表しています．この数理モデルは最も簡単な捕食者・非捕食者モデルとなっています．このとき，X 種，N種の初期個体群密度 x(0), n(0) をそれぞれ x0 , n0 とすると，(10) の 2 式からd(x + βn) = 0dt

を得ることができ

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

Y + N −→ (1 + βy )Yという式が得られます． このとき数理モデルはdx= αx βx xn,dt(12)dy= αy βy yn,dt

(13)

dn
= −αx xn − αy yndt

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

dx

= ax (1 − bx x)x − kx xy,
dt dy = a (1 − b y)y − k xyyyydt(16)となります．ただし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 )βxβyαy, by =, ky = αy .
βx
βy x0 + βx y0 + βx βy n0βx

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

3. まとめ
本稿では，数理モデリングを行うための基本原理を説明し，その原理に基づいて非常にシンプルな仮定での数理モデリングの例を説明し，化学反応モデル，個体群動態モデルの構成について説明しました．本稿では説明しませんでしたが，個体群動態のほかに感染症モデルも化学反応モデルと似た形式で記述できます [9]．また，生体内で生じる生化学反応（例えば，酵素分解反応）や細胞内でのカルシウムイオン興奮なども同じような方法で導出することができます [10]．市民講演会では，ここで説明した数理モデル以外に，ろうそく火炎振動子の同期現象 [11] やメトロノームの同期現象に対する数理モデリングについて実験動画とシミュレーション動画を交えて紹介しました．
そして，最近研究している表皮の数理モデリングについても紹介しました [12] が，これらの説明はここでは割愛させていただきます．最後になりますが，2 つのことを強調しておきたいと思います．1 つは，与えられた微分方程式の変数や定数に現象の言葉として意味を持てば，その微分方程式は数理モデルになるということです．微分方程式の研究が何の役に立つの？ とよく聞かれますが，研究対象とする微分方程式が現象としての意味を持てばその答えは明白だと思います．もう 1 つは，数理モデリングとは数式による現象の模倣であるということです．基本原理がわからない現象に対して，基本原理の一端を理論的に捉えるための道具なのです．基本的に仮説を与える時点で主観が入るため，どうしても数理モデリングは客観性を失います．従って，数理モデリングによって現象を再現できたからといって現象の数理モデルができたとは言えないのです．得た仮説から可能な限り主観を取り除く努力をすることが必要なので，数理モデルから示唆されたメカニズムや予測を実験で検証することが重要になってきます．
市民講演会での講演を聴いて，一人でも現象の数理モデリングに興味を持っていただ
けたなら幸いです．

参考文献[1] “微分方程式で数理モデルをつくろう”，デヴィッド・バージェス，モラル・ボリー，日本
評論社[2] “非平衡系の科学 III”，三池 秀敏，森 義仁，山口 智彦，講談社.
[3] B. P. Belousov, in Collections of abstracts on radiation medicine(Medgiz, Moscow), 1958,145–.
[4] A. M. Zhabotinsky, Biofizika, 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