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

​