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 のみであることが分かります。
そのためプログラム中のコメントに示したような処理を行いました。