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")

​