R で数理生物学入門：種間競争のダイナミクスとアイソクライン法

今回から第 2 章に入ります。
数理生物学入門 p.13-17 のロトカ・ボルテラ競争系です。
これは以下の式

$dx/dt=r_1x(1-(x+ay)/K_1)$
$dy/dt=r_2y(1-(bx+y)/K_2)$

で表されます。
この式において、$a$ は $y$ 一個体分の影響を $x$ 一個体分の影響に換算する係数で、$b$ はその逆になっています。$ay, bx$ がなければロジスティック成長の式と同じです。つまり、ロジスティック成長を基に、相手の個体数を考慮した式になっているということです。詳しい説明については本書をご参照ください。

ここではプログラムについて説明します。
以下がロトカ・ボルテラ競争系のプログラムです。

#作成日時 : 13:55 2016/12/03
#作成者名：中野　善
rm(list=ls())
gc();gc()

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

# 変数の設定 ---------------------------
x = 1            # 個体数の初期値
y = 1
dt = 0.05         # ちょっとだけ進める時間
T = 1000          # 計算回数

r1=1
r2=2
K1=80
K2=70
a=0.7
b=1.1

# 計算 --------------------------------
Nx <- x
Ny <- y

for(i in 1:T){
dx <-  r1*x*(1-(x+a*y)/K1)*dt
dy <-  r2*y*(1-(b*x+y)/K2)*dt
x <- x+dx
y <- y+dy
Nx <- c(Nx, x)
Ny <- c(Ny, y)
}

plot(Nx, Ny, type="l", xlim=c(0, 100), ylim=c(0, 100))
segments(0, K2, K2/b, 0, lty=2)
segments(0, K1/a, K1, 0, lty=2)

​

式が 2 つになって、ややこしくなったような気がするのですが、よく見るとこれまでの式の拡張で、それほどややこしくないんです。 for のすぐ下の 2 行はロトカ・ボルテラ競争系の式 $dt$ を右辺に移項して表示しただけです。これはロジスティック成長の式に $ay, bx$ を追加したものになってます。そして、次の 2 行で x, y を更新して、さらに次の 2 行で Nx, Ny に計算結果を入れています。このプログラムの最初の「変数の設定」の x, y を変更することで、個体数の初期値を変更できるようになっています。



上記プログラムの個体数の初期値を複数準備し、繰返し計算させたものが次のプログラムです。このプログラムで図 2.1 を再現できます。

#作成日時 : 13:55 2016/12/03
#作成者名：中野　善
rm(list=ls())
gc();gc()

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

# 変数の設定 ---------------------------
intx = c(1, 0.1, 5, 10, 1, 50, 90, 100)            # 個体数の初期値
inty = c(40, 0.1, 2, 1, 100, 100, 100, 50)
dt = 0.05         # ちょっとだけ進める時間
T = 1000          # 計算回数

r1=1
r2=2
K1=80
K2=70
a=0.7
b=1.1

# 計算 --------------------------------
NX <- NY <- as.data.frame(matrix(0, nrow=T+1, ncol=length(intx)))

for(j in 1:length(intx)){
Nx <- x <- intx[j]
Ny <- y <- inty[j]

for(i in 1:T){
dx <-  r1*x*(1-(x+a*y)/K1)*dt
dy <-  r2*y*(1-(b*x+y)/K2)*dt
x <- x+dx
y <- y+dy
Nx <- c(Nx, x)
Ny <- c(Ny, y)
}    # for(i in 1:T)

NX[,j] <- Nx
NY[,j] <- Ny
}    # for(j in 1:length(intx))

i=1
plot(NX[,i], NY[,i], type="l", xlim=c(0, 100), ylim=c(0, 100), xlab="第 1 種 x", ylab="第 2 種 y", lab=c(3,3,0), bty="l")
segments(0, K2, K2/b, 0, lty=2)
segments(0, K1/a, K1, 0, lty=2)

for(i in 1:ncol(NY)){
lines(NX[,i], NY[,i])
}

points(K1, 0, pch=16)
points(0, K2, pch=16, col="white"); points(0, K2)

text(K2/b, 0, label="K2/b", pos=2, cex=0.8)
text(K1, 0, label="K1", pos=2, cex=0.8)

text(0, K2, label="K2", pos=4, cex=0.8)
text(0, K1/a, label="K1/a", pos=4, cex=0.8)

legend("topright", legend="(a)", bty="n")

#作成日時 : 13:55 2016/12/03
#作成者名：中野　善
rm(list=ls())
gc();gc()

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

# 変数の設定 ---------------------------
intx = c(1, 2, 10, 30, 25, 70, 100, 100, 100)            # 個体数の初期値
inty = c(40, 2, 2, 2, 100, 100, 100, 50, 25)
dt = 0.05         # ちょっとだけ進める時間
T = 1000          # 計算回数

r1=1
r2=2
K1=75
K2=80
a=0.7
b=0.8

# 計算 --------------------------------
NX <- NY <- as.data.frame(matrix(0, nrow=T+1, ncol=length(intx)))

for(j in 1:length(intx)){
Nx <- x <- intx[j]
Ny <- y <- inty[j]

for(i in 1:T){
dx <-  r1*x*(1-(x+a*y)/K1)*dt
dy <-  r2*y*(1-(b*x+y)/K2)*dt
x <- x+dx
y <- y+dy
Nx <- c(Nx, x)
Ny <- c(Ny, y)
}    # for(i in 1:T)

NX[,j] <- Nx
NY[,j] <- Ny
}    # for(j in 1:length(intx))

i=1
plot(NX[,i], NY[,i], type="l", xlim=c(0, 100), ylim=c(0, 100), xlab="第 1 種 x", ylab="第 2 種 y", lab=c(3,3,0), bty="l")
segments(0, K2, K2/b, 0, lty=2)
segments(0, K1/a, K1, 0, lty=2)

for(i in 1:ncol(NY)){
lines(NX[,i], NY[,i])
}

points(K1, 0, pch=16, col="white"); points(K1, 0)
points(0, K2, pch=16, col="white"); points(0, K2)

# maxima: solve([y=-K2/(K2/b)*x+K2, y=(-K1/a)/K1*x+K1/a], [x,y]);
points((a*K2-K1)/(a*b-1), -(K2-b*K1)/(a*b-1), pch=16)

text(K2/b, 0, label="K2/b", pos=2, cex=0.8)
text(K1, 0, label="K1", pos=2, cex=0.8)

text(0, K2, label="K2", pos=4, cex=0.8)
text(0, K1/a, label="K1/a", pos=4, cex=0.8)

legend("topright", legend="(b)", bty="n")

#作成日時 : 13:55 2016/12/03
#作成者名：中野　善
rm(list=ls())
gc();gc()

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

# 変数の設定 ---------------------------
intx = c(0.5, 10, 12, 15, 25, 40, 60, 100, 100, 100)            # 個体数の初期値
inty = c(0.1, 5, 5, 2, 100, 100, 100, 100, 50, 20)
dt = 0.05         # ちょっとだけ進める時間
T = 1000          # 計算回数

r1=1
r2=2
K1=90
K2=90
a=1.4
b=1.5

# 計算 --------------------------------
NX <- NY <- as.data.frame(matrix(0, nrow=T+1, ncol=length(intx)))

for(j in 1:length(intx)){
Nx <- x <- intx[j]
Ny <- y <- inty[j]

for(i in 1:T){
dx <-  r1*x*(1-(x+a*y)/K1)*dt
dy <-  r2*y*(1-(b*x+y)/K2)*dt
x <- x+dx
y <- y+dy
Nx <- c(Nx, x)
Ny <- c(Ny, y)
}    # for(i in 1:T)

NX[,j] <- Nx
NY[,j] <- Ny
}    # for(j in 1:length(intx))

i=1
plot(NX[,i], NY[,i], type="l", xlim=c(0, 100), ylim=c(0, 100), xlab="第 1 種 x", ylab="第 2 種 y", lab=c(3,3,0), bty="l")
segments(0, K2, K2/b, 0, lty=2)
segments(0, K1/a, K1, 0, lty=2)

for(i in 1:ncol(NY)){
lines(NX[,i], NY[,i])
}

points(K1, 0, pch=16)
points(0, K2, pch=16)

# maxima: solve([y=-K2/(K2/b)*x+K2, y=(-K1/a)/K1*x+K1/a], [x,y]);
points((a*K2-K1)/(a*b-1), -(K2-b*K1)/(a*b-1), pch=16, col="white")
points((a*K2-K1)/(a*b-1), -(K2-b*K1)/(a*b-1))

text(K2/b, 0, label="K2/b", pos=2, cex=0.8)
text(K1, 0, label="K1", pos=2, cex=0.8)

text(0, K2, label="K2", pos=4, cex=0.8)
text(0, K1/a, label="K1/a", pos=4, cex=0.8)

legend("topright", legend="(c)", bty="n")

​

R で数理生物学入門：スイッチする捕食者

今回は数理生物学入門 p.9-10 のスイッチする捕食者についてです。
これは捕食の影響が餌生物の個体数によって変わってくるという現象を表したもので、次の数式で表されます。

$dx/dt=rx(1-x/K)-Px^2/(1+x^2)$　...(1)

第一項がロジスティック式で増殖を表していて、第二項が捕食による減少量を表しています。この各項をグラフにしたものが図 1.6 で、以下のプログラムで再現できます。

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

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

# 変数の設定 ---------------------------
K = 20
r = 0.5
P = 1.15

x = seq(0,K+1, by=0.1)
dx.dt1 <- r*x*(1-x/K)            # dx/dt の第 1 項：増殖
dx.dt2 <- P*(x^2)/(1+x^2)        # dx/dt の第 2 項：捕食

plot(x, dx.dt1, type="l", ylab="害虫個体数の時間変化 dx/dt", xlab="餌生物の個体数 x")
lines(x, dx.dt2)
abline(h=0)

# 平衡点 A, B, C の値を求める --------------
act1 <- dx.dt1-dx.dt2
act2 <- (act1+0.001)/abs(act1+0.001)
act3 <- x[which(diff(act2)!=0)+1]

A <- act3[1]
B <- act3[2]
C <- act3[3]

points(c(0, A, B, C), c(0,0,0,0), pch=16, col=c("white",1,"white",1))
points(c(0, A, B, C), c(0,0,0,0), pch=c(21,16,21,16))

text(c(0, A, B, C), c(0,0,0,0), label=c("0", "A", "B", "C"), pos=1, cex=0.8)
text(K, max(dx.dt1)/5, label="増殖")
text(K, max(dx.dt2), label="捕食", pos=3)

​

図 1.7 については以下のプログラムを実行すると再現できます。図 1.7 は図 1.6 の増殖と捕食の交点 A, B, C の値が、捕食の強さ $P/r$ に依存して変化する挙動を表しているようです。
すなわち

$rx(1-x/K)-Px^2/(1+x^2)=0$

となる $x$ を求めてやれば良いわけです。Maxima を使って解を求めたところ、自分のやり方がまずかったのか、虚数が出てしまいました。これでは図示するのが大変そうなので、これも数値的に解くことにしました。それが以下のプログラムです。

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

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

# 変数の設定 ---------------------------
K = 20
r = 0.5
P <- seq(0.01, 4, length=1000)

A <- B <- C <- NULL
x <- seq(0,K+1, by=0.001)

for(i in 1:length(P)){
dx.dt1 <- r*x*(1-x/K)                # dx/dt の第 1 項：増殖
dx.dt2 <- P[i]*(x^2)/(1+x^2)         # dx/dt の第 2 項：捕食

act1 <- dx.dt1-dx.dt2
act2 <- (act1+0.001)/abs(act1+0.001)
act3 <- x[which(diff(act2)!=0)+1]

A <- c(A, act3[1])
B <- c(B, act3[2])
C <- c(C, act3[3])

}    # for(i in 1:length(P))

# -------------------------------------------------------------------------
# 解が一つしかなく、その値が大きい場合、その解は C である。
# しかし上記計算ではその解は A に入ってしまうため、ここで C に移している。

C[which(A > K/2)] <- A[which(A > K/2)]
A[which(A > K/2)] <- NA

# --------------------------------------------------------------------------

plot(P/r, A, type="l", ylim=c(0,K), ylab="平衡状態での個体数 x", xlab="捕食圧の強さ P/r")
lines(P/r, B, lty=2)
lines(P/r, C)

text(mean(P[which(!is.na(A))])/r, mean(A, na.rm=T), label="安定 A", pos=3)
text(mean(P[which(!is.na(B))])/r, mean(B, na.rm=T), label="不安定 B", pos=2)
text(mean(P[which(!is.na(C))])/r, mean(C, na.rm=T), label="安定 C", pos=2)

​

上記プログラムでは $x$ を 0 から $K+1$ まで 0.001 刻みの値を設定しています。ここで、計算した act1 が 0 になったときの $x$ を取ってこようと考えたのですが、都合よく 0 になってくれるとも限らないので、act1 の正負の符号が変わったときの $x$ を取ってくることにしました。

そこで変数 act2 で act1 を -1 と 1 に変換し（+0.001 は act1 が 0 のときに NaN が出てくるのを防ぐためのもの）、変数 act3 で差分を取ることで、符号が変わったところを 2 または -2 に変換しました（符号が同じところは 0）。そして 2 または -2 のところの要素番号を取ってきて（差分を取ったので要素番号に 1 を足して）、その要素番号に対応する $x$ の値を抽出しました。そしてその値を A, B, C に振り分けました。

これを 0.01 から 4 まで 1000 分割した $P$ を使って for で反復計算し、それぞれの $P$ に対応する A, B, C の値をまとめました。



上記プログラム中のコメント「解が一つしかなく、その値が大きい場合、その解は C である。」について捕捉のプログラムを以下に示します。

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

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

K = 20
r = 0.5
P = c(0.01, 0.92, 1.8, 2.55, 4)

i=1
x = seq(0,K+1, by=0.1)
dx.dt1 <- r*x*(1-x/K)                # dx/dt の第 1 項：増殖
dx.dt2 <- P[i]*(x^2)/(1+x^2)         # dx/dt の第 2 項：捕食

act1 <- dx.dt1-dx.dt2
plot(x, act1, type="l", ylim=c(-4, 5)); abline(h=0)
#plot(x, act1, type="l", ylim=c(-0.1, 0.1), xlim=c(0,0.5)); abline(h=0)

for(i in 1:length(P)){
x = seq(0,K+1, by=0.1)
dx.dt1 <- r*x*(1-x/K)
dx.dt2 <- P[i]*(x^2)/(1+x^2)

act1 <- dx.dt1-dx.dt2
lines(x, act1)
}

​

この図では P の値によって解（act1=0）の数が変わることが見て取れます。
一番上の曲線では解は 0 以外では 1 つであり、曲線が下がってくる（すなわち $P$ が大きくなる）にしたがって解の数が 2 つ（0 に接している）、3 つとなっています。さらに曲線が下がってくると、解の数は再び 2 つになり、やがて 1 つになります。
このとき $P$ が小さいとき（一番上の曲線）の解は C のみであり、$P$ が大きいとき（一番下の曲線）の解は A のみであることが分かります。
そのためプログラム中のコメントに示したような処理を行いました。

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)

​

R で数理生物学入門：ロジスティック成長

今回は数理生物学入門 p.4 のロジスティック成長です。

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

$r$ は内的自然増加率、$K$ は環境収容力を表しています。
環境収容力は、その環境に何個体入れるかという値です。
つまり、個体数 $x$ が既に $K$ に達していると、それ以上増えることができないということになります。
実際に $x=K$ とすると、右辺は $rK(1-K/K)=0$ となり、すなわち $dx=0$ となって、確かに増加しないことがわかります。
逆に $x$ の値が小さいときは $r$ の値がほぼダイレクトに反映されるようになっていることがわかります。

数式の解説は本書に任せるとして、この場ではプログラムを考えてみます。
早速ですが、この数式のプログラムは以下の通りです。

#作成日時 : 15:33 2016/10/15
#作成者名：T. Nakano
rm(list=ls())
gc();gc()

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

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

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

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

​

このプログラムについて説明します。
と言っても、基本的には前回の指数増殖と同じで、

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

に持ち込んでいます。
ここで $dx$ は式(1)から、

$dx=rx(1-x/K)\times dt$　…(3)

となります。
今回の式 (3) は前回の式 (3) と比較するとわかるように、マルサス係数 $m$ を $r(1-x/K)$ に置き換えたものになっているので、プログラムの方もそのように変更しました。




以下のプログラムで本書のグラフを再現できます。

#作成日時 : 15:33 2016/10/15
#作成者名：T. Nakano
rm(list=ls())
gc();gc()

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

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

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

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

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

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

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

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

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

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

text(10, 100, label="x = K", pos=3)

​

shiny を使うと以下のようになります。

ui.R

library(shiny)

shinyUI(fluidPage(

  # Application title
  titlePanel("ロジスティック成長"),

  sidebarLayout(
    sidebarPanel(
      numericInput("r","r =", value = 1),
      numericInput("K","K =", value = 100),
      numericInput("x","x =", value = 1),
      numericInput("dt","dt =", value = 0.1),
      numericInput("T","T =", value = 100)
    ),

    mainPanel("ロジスティック成長",
      plotOutput("logi")
    )
  )
))

​

server.R

library(shiny)

shinyServer(function(input, output) {

  output$logi <- renderPlot({
    # 変数の設定 ---------------------------
    r = input$r            # 内的自然増加率
    K = input$K            # 環境収容力
    x = input$x            # 個体数の初期値
    dt = input$dt          # ちょっとだけ進める時間
    T = input$T            # 計算回数

    ymax <- 150
    if(ymax < x) ymax <- x

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

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

  })

})

​