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)

  })

})

​

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

​

番外編：お菓子300円問題

突然ですが…、
小学校の遠足の時「お菓子は 300 円までね～」っていうルールは無かったですか？そして、これをキッチリ使い切りたいっていう悩ましい問題に陥って駄菓子屋さんの前で長い時間を過ごしたことかと思います。ということで、今回は小学生の気持ちに戻って、「お菓子 300 円問題」を R で解決してみたいと思います。

計算自体はとても簡単です。
予め、欲しいお菓子の名前と金額を csv ファイルにまとめておき、ランダムに取り出してお小遣いから引いていきます。そして、お小遣いが0円になったら計算を終了します。もしお小遣いがマイナスになったら初めから計算し直します。


使い方 -----------------
・欲しいお菓子の名前と金額を csv ファイルを作成します。例はこちら。
・以下のコードを実行し、csv ファイルを選択します。
以上

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

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

# 変数の設定 -----------------------------------------------------
# お小遣い
bgt <- 300    # 円

# データの読込み -------------------------------------------------
library(tcltk)
fileName <- tclvalue(tkgetOpenFile())
act1 <- unlist(strsplit(fileName,"/"))
act1 <- act1[-length(act1)]

directory <- NULL
if(length(act1)!=1){
for(i in 1:(length(act1)-1)){
directory <- paste(directory,act1[i],"/",sep="")
}
}
directory <- paste(directory,act1[length(act1)],sep="")
setwd(directory)
DT <- read.csv(fileName, header=T)

# 計算開始 ------------------------------------------------------
N <- round(bgt/min(DT$金額))

for(j in 1:1000){
WR <- NULL
for(i in 1:N){
elm <- round(1+runif(1)*(nrow(DT)-1))
WR <- rbind(WR, DT[elm,])
if(bgt-sum(WR$金額)<=0) break
}
if(bgt==sum(WR$金額)){
WR <- WR[order(WR$お菓子),]
write.csv(WR, file="お菓子購入リスト.csv", row.names=F)
View(WR)
break
}
}

# コメント --------------------
if(j==1000){
cat(paste(bgt,"円ちょうどになりませんでした。", "\n","「お菓子と金額.csv」を変更して再計算してください。", sep=""), "\n")
}else if(j < 1000){
cat(paste(bgt,"円ちょうどになりました。", "\n", "ファイルを保存しましたのでお菓子を買いに行きましょう。", sep=""), "\n")
}

​

今回は「お菓子300円問題」ということでしたが、例えば、友達とバーベキュー等をするときにも使えるかと思います。バーベキュー等のときは駄菓子と違って商品の金額のキリが悪いことがあります（牛肉876円とか）。そういうときは一の位が 1 円か 9 円のものが入っている方が計算がうまくいきますよ。