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)

  })

})

​