R で数理生物学入門：指数増殖

巌佐庸先生の数理生物学入門は非常に面白く、分かりやすく書かれており、感銘を受けました。この素晴らしい書籍を手にされた読者のうち、数式が苦手な方の一助になればと思い、自分が理解を深めるために作った R コードを公開したいと思いました。

というわけで、本書で紹介されている数式を R に起こして勝手にサポートしたいと思います！具体的には、数値解を求める R コードを公開します (どこまでできるか分かりませんが)。解析解 (数式としての解) の説明については、本書や他の書籍をご参照下さい。

まずは手始めに本書 p.3 の

$dx/dt = mx$　　　…　(1)

についてです。
本書によると

$x(t) = x(0)e^{mt}$ が解であることは式の両辺に直接代入して確かめることができる。
この結果は、一定の環境がしばらく続くと個体数が時間とともに指数関数を描いて増大することを表している。

とのことですが、数式が苦手だと、いきなり「？？？」が浮かぶことになります。数式が得意な人達には分からない感覚だと思いますが、現に「？？？」が浮かびフリーズすることがあるんです。

でも、諦めなくていいんです！
R を使えば、別の道から解いていけるんです。

本当なら積分して解析解を求めるのがいいんだと思うんですが、数学的テクニックがなければそれも叶いません。そういうときは数値解を計算していくことで近似的に解を求めることができるんです。ここでいう数値解を計算するということは、具体的な値を代入して、数値として解を算出していくことです。今回の場合だと、式 (1) は $x$ の増加を記述するための数式なので、数値解では具体的な $x$ の値を算出していくことが最終目標になります。

今回が基本の数式となるので、ちょっと真面目に考えてみます。
説明が少しまどろっこしくなりますが、ご容赦願います。

式 (1) の左辺 $dx/dt$ は、言い換えると「$t$ がちょっと動いたときに $x$ がどれぐらい動いたか」を表しています。「$t$ がちょっと動いたとき」が $dt$ で、「$x$ がどれぐらい動いたか」が $dx$ で表されているわけです。もう少し格好良く言えば、それぞれ $t$ の変化量、$x$ の変化量と表現できます。

本書より、$x$ が個体数、$t$ が時刻なので、式 (1) の左辺 $dx/dt$ は、言い換えると「時間がちょっと進んだときの個体数の増加分 ($m < 0$ なら減少分) 」を表していることになります。「時間がちょっと進んだ」というのは、単位を任意で 1 進んだとすると、$t+1$ と表せます。$t$ を今の時刻とすると、$t+1$ は次の時刻と言えるわけです。

ここで、「次の時刻の個体数」というものは、「今の時刻の個体数」に「時間がちょっと進んだときの個体数の増加分」を加えたものになります。なので、次の時刻の個体数を $x_{t+1}$、今の時刻の個体数を $x_t$、時間がちょっと進んだときの個体数の増加分を $dx$ とすると、次の時刻の個体数は、

$x_{t+1}=x_t+dx$　　… 　(2)

と表せます。そして $dx$ は式 (1) を変形すると、

$dx=mx\times dt$　　…　(3)

というように求めることができます。ここで式 (3) の $x$ は「今の個体数」になるので $x_t$ となります。ここまでくれば、式 (2) と (3) から $x_{t+1}$ を求めることができます。そして求めた $x_{t+1}$ を $x_t$ に代入し、次の $x_{t+1}$ を順々に求めていけばいいわけです。

というわけで、この方針でプログラムを組んだものが以下のコードです。

#作成日時 : 13:14 2016/09/25
#作成者名：T. Nakano
rm(list=ls())
gc();gc()

# /////////////////////////////////////////////////////////////////////

# 変数の設定 ---------------------------
m = 0.1            # マルサス係数
x = 10             # 個体数の初期値
dt = 0.1           # ちょっとだけ進める時間

# 計算 --------------------------------
N <- x
for(i in 1:1000){
dx <-  m*x*dt      # 式(3)
x <- x+dx          # 式(2)
N <- c(N, x)
}

plot(0:1000*dt, N, type="l", xlab="t")

​