R で数理生物学入門：低密度の影響

今回は数理生物学入門 p.7 の低密度の影響についてです。
いわゆるアリー効果というやつで、次の式で表されます。

$dx/dt = rx(x-a)(1-x/K)$　... (1)

$a$ は低密度の影響が現れる閾値で、$0<a<K$ と定義されます。$x$ がこの閾値を超えるか超えないかで個体数の挙動が変わってくるというわけです。試しに $x<a$ を考えてみると、式 (1) の $(x-a)$ がマイナスになるので $dx/dt$ もマイナスとなります。変化量がマイナスなので確かに個体数は時間とともに減少します。逆に $x>a$ を考えてみると、$dx/dt > 0$ で個体数が増加することがわかります。



今回も詳しい解説は本書を参照してもらうとして、ここでは式 (1) を表現するプログラムを示します。以下がそのプログラムです。

#作成日時 : 11:07 2016/11/05
#作成者名：T. Nakano
rm(list=ls())
gc();gc()

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

# 変数の設定 ---------------------------
r = 0.03            # 内的自然増加率
K = 100             # 環境収容力
x = 30              # 個体数の初期値
a = 20              # 0 < a < K
dt = 0.1            # ちょっとだけ進める時間
T = 100             # 計算回数

# 計算 --------------------------------
N <- x
for(i in 1:T){
dx <-  r*x*(x-a)*(1-x/K)*dt
x <- x+dx
N <- c(N, x)
}

# 描画 -------------------------------
plot(0:T*dt, N, type="l", xlab="t")
abline(h=K, lty=2)

​

このプログラムはロジスティック成長のプログラムとほとんど同じです。ロジスティック成長のときにマルサス係数 $m$ を $r(1-x/K)$ に置き換えたのですが、今回はこれに $(x-a)$ を書き加えただけです。そして、$a$ が入ったので変数の設定で $a$ を定義しておきました。




以下のプログラムで本書の図1.5(a) を再現できます。

#作成日時 : 11:07 2016/11/05
#作成者名：T. Nakano
rm(list=ls())
gc();gc()

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

# 変数の設定 ---------------------------
r = 0.03            # 内的自然増加率
K = 100             # 環境収容力
x = 19              # 個体数の初期値
a = 20              # 0 < a < K
dt = 0.1            # ちょっとだけ進める時間
T = 100             # 計算回数

# 計算 --------------------------------
N <- x
for(i in 1:T){
dx <-  r*x*(x-a)*(1-x/K)*dt
x <- x+dx
N <- c(N, x)
}

plot(0:T*dt, N, type="l", xlab="時間 t", ylab="個体数 x", ylim=c(0, 160))

# 変数の設定 ---------------------------
x = 21             # 個体数の初期値

# 計算 --------------------------------
N <- x
for(i in 1:T){
dx <-  r*x*(x-a)*(1-x/K)*dt
x <- x+dx
N <- c(N, x)
}

lines(0:T*dt, N, type="l")

# 変数の設定 ---------------------------
x = 150             # 個体数の初期値

# 計算 --------------------------------
N <- x
for(i in 1:T){
dx <-  r*x*(x-a)*(1-x/K)*dt
x <- x+dx
N <- c(N, x)
}

lines(0:T*dt, N, type="l")

abline(h=c(K, a), lty=2)
text(10, K, label="x = K", pos=3)
text(10, a, label="x = a", pos=3)

​

また、以下のプログラムで本書の図1.5(b) を再現できます。これは個体数と増殖率の関係なので、$rx(x-a)(1-x/K)$ の $x$ に 0 から 120 までの値を与えて計算した結果をグラフにしました。

#作成日時 : 11:07 2016/11/05
#作成者名：T. Nakano
rm(list=ls())
gc();gc()

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

# 変数の設定 ---------------------------
r = 0.03            # 内的自然増加率
K = 100             # 環境収容力
a = 20              # 0 < a < K

# 計算 --------------------------------
x = 0:120
dx.dt <-  r*x*(x-a)*(1-x/K)

plot(x, dx.dt, type="l", ylim=c(-40, 40), ylab="増殖率", xlab="個体数 x")
points(c(0, a, K), c(0, 0, 0), pch=c(16, 21, 16))
abline(h=0)

text(a, 0, label="a", pos=1)
text(K, 0, label="K", pos=1)

​