水産資源解析学

Lecture 06: 種内および種間競争について

Greg Nishihara (長崎大学)

2012 November 15

ここに記載したRコードを実行するためには、次のパッケージのインストールが必要です。パッケージをインストールする場合は「install.packages(“packagename”)」を実行してください。例えば、deSolveというパッケージをインストールする場合には、「install.packages(“deSolve”)」を実行しましょう。

require(deSolve)

## Loading required package: deSolve

require(ggplot2)

## Loading required package: ggplot2

require(plyr)

## Loading required package: plyr

require(reshape2)

## Loading required package: reshape2

幾何級数的個体群成長 (離散世代)

\( \lambda \) は単位時間あたり増殖率、\( b \) と \( d \) 出生率と死亡率、\( N \) は個体数。 \[ \begin{aligned} N_{t+1} &= \lambda N_t \\ \lambda &= b-d \\ N_{t+1} &= (b-d) N_t \\ \end{aligned} \] では、まずはgeometric_growth() という、関数を定義する。

geometric_growth = function(N, b, d) {
    (b - d) * N
}

初期パラメータは、 \( b=5.3 \) と \( d=5 \)、初期個体数は \( N_0=10 \)。

N0 = 10  # 初期個体数
b = 5.3  # 出生率
d = 5  # 死亡率
Time = 0:5  # 単位時間
N = rep(N0, length(Time))
for (i in 1:(length(N) - 1)) {
    N[i + 1] = N[i] + geometric_growth(N[i], b = b, d = d)
}

Parameters

\( N_0= \) 10; \( b= \) 5.3; \( d= \) 5

作図にはggplot2というライブラリーを活用するが、まずはdata.frameを構築します。

dset = data.frame(Time = Time, N = N)
dset$lambda = "> 1"

dset の先頭の値を取り出して表示する。

head(dset)

##   Time     N lambda
## 1    0 10.00    > 1
## 2    1 13.00    > 1
## 3    2 16.90    > 1
## 4    3 21.97    > 1
## 5    4 28.56    > 1
## 6    5 37.13    > 1

ggplot()の関数を使って、作図。

ggplot(data = dset, aes(x = Time, y = N, color = lambda)) + geom_line(size = 4, 
    alpha = 0.6) + labs(x = "Time", y = "Population size") + labs(color = expression(lambda)) + 
    theme_bw(20)

plot of chunk unnamed-chunk-6

\( \lambda \) を変えるとこのようになる。

N0 = 10
b2 = 5.3
b1 = 5
b0 = 4.5
d = 5
Time = 0:5
N2 = rep(N0, length(Time))
N1 = rep(N0, length(Time))
N0 = rep(N0, length(Time))
for (i in 1:(length(N0) - 1)) {
    N2[i + 1] = N2[i] + geometric_growth(N2[i], b = b2, d = d)
    N1[i + 1] = N1[i] + geometric_growth(N1[i], b = b1, d = d)
    N0[i + 1] = N0[i] + geometric_growth(N0[i], b = b0, d = d)
}

dset = data.frame(Time = rep(Time, 3), N = c(N2, N1, N0), lambda = factor(c(rep(1, 
    length(N2)), rep(2, length(N1)), rep(3, length(N0))), label = c("> 1", "= 1", 
    "< 1")))

ggplot(data = dset, aes(x = Time, y = N, color = lambda)) + geom_line(alpha = 0.6, 
    size = 2) + labs(x = "Time", y = "Population size", color = expression(lambda)) + 
    theme_bw(20)

plot of chunk unnamed-chunk-7

指数的個体群成長 (連続世代)

\( r \) は瞬間増殖率、\( b \) と \( d \) 出生率と死亡率、\( N \) は個体数。 \[ \begin{aligned} \frac{dN}{dt} &= r N \\ r &= b-d \\ \frac{dN}{dt} &= (b-d) N \end{aligned} \] ではexponential_growth() という、関数を定義しましょう。この場合、deSolveの連立微分方程式のソルバーを使うので、関数の書き方を工夫します。

exponential_growth = function(Time, State, Param) {
    with(as.list(c(State, Param)), {
        dN = (b - d) * N
        return(list(c(dN)))
    })
}

ソルバーに必要な初期値はここで定義します。初期個体数は \( N_0=10 \)、出生率と死亡率は \( b=5.3 \) と \( d=5 \) です。

State = c(N = 10)
Param = c(b = 5.3, d = 5)
Time = seq(0, 5, by = 0.001)
model = ode(State, Time, exponential_growth, Param, method = "ode45")
dataout = as.data.frame(model)

ここでは、reshape2のパッケージにあるmelt()という関数をつかって、dataoutをもっと都合にいいフォーマットにします。

dset = melt(dataout, id = "time", variable.name = "Species")
dset$r = factor(1, label = "> 0")

Parameters

\( N_0= \) 10; \( b= \) 5.3; \( d= \) 5

結果

ggplot(data = dset, aes(x = time, y = value, color = r)) + geom_line(alpha = 0.6, 
    size = 2) + labs(x = "Time", y = "Population size") + theme_bw(20)

plot of chunk unnamed-chunk-11

\( r \) を帰るとこのような結果になる。

State = c(N = 10)
Time = seq(0, 5, by = 0.001)

Param2 = c(b = 5.3, d = 5)
Param1 = c(b = 5, d = 5)
Param0 = c(b = 4.5, d = 5)
model2 = ode(State, Time, exponential_growth, Param2, method = "ode45")
model1 = ode(State, Time, exponential_growth, Param1, method = "ode45")
model0 = ode(State, Time, exponential_growth, Param0, method = "ode45")
dataout2 = as.data.frame(model2)
dataout1 = as.data.frame(model1)
dataout0 = as.data.frame(model0)
dset2 = melt(dataout2, id = "time", variable.name = "r")
dset1 = melt(dataout1, id = "time", variable.name = "r")
dset0 = melt(dataout0, id = "time", variable.name = "r")
dset2$r = factor(1, label = "> 0")
dset1$r = factor(2, label = "= 0")
dset0$r = factor(3, label = "< 0")
dset = rbind(dset2, dset1, dset0)
ggplot(data = dset, aes(x = time, y = value, color = r)) + geom_line(alpha = 0.6, 
    size = 2) + labs(x = "Time", y = "Population size") + theme_bw(20)

plot of chunk unnamed-chunk-12

Lotka-Volterra Model 1

被食者と捕食者の相互作用

次は、被食者と捕食者の相互作用を成長モデルに加えます。被食者の個体数 (N) の変化は捕食者の個体数 (P) と結ばれていて、同じように捕食者の個体数は被食者の個体数と結ばれている。 \[ \begin{aligned} \frac{dP}{dt} &= f(N, P) \\ \frac{dN}{dt} &= f(P, N) \end{aligned} \] まずは、捕食者の出生率は指数関数的に増加させる。死亡率は捕食圧に依存する。出生率は \( r_N \times N \) 、捕食圧は \( \alpha_{NP}\times N\times P \) です。\( \alpha_{NP} \) は捕食圧の効率を表している。 \[ \frac{dN}{dt}=r_N N - \alpha_{NP} N P \] 捕食者の出生率は餌を得る効率と餌を個体数に変換する効率と結びつける。よって、 \( \alpha_{PN} \alpha{NP} \times P \ times N \) となる。\( \alpha_{PN} \) は餌の変換効率である。死亡率は捕食者の個体数の一定の割合 \( r_P \times P \) で増加させる。 \[ \frac{dP}{dt}=\alpha_{PN} \alpha_{NP} P N - r_P P \]

この連立微分方程式はdeSolveのode()関数を使って結果を求めます。

LV1 = function(Time, State, Param) {
    with(as.list(c(State, Param)), {
        dN = rN * N - alphaNP * N * P
        dP = alphaNP * alphaPN * N * P - rP * P
        return(list(c(dN, dP)))  # return()の定数の順位は必ず連立微分方程式の順位と同じ順位にする。
    })
}

例１。

State = c(N = 10, P = 20)
Param = c(rP = 1, rN = 2, alphaNP = 0.8, alphaPN = 0.5)
Time = seq(0, 50, by = 0.01)
model = ode(State, Time, LV1, Param, method = "ode45")
dataout = as.data.frame(model)

dset = melt(dataout, id="time", variable.name="Species")
dset$Species = as.factor(dset$Species)

モデルのパラメーター

\( N_{0}= \) 10; \( P_{0}= \) 20; \( r_P= \) 1; \( r_N= \) 2; \( \alpha_{PN}= \) 0.5; \( \alpha_{NP}= \) 0.8

データの先頭と末尾を取り出して表示。

rbind(head(dset), tail(dset))

##        time Species   value
## 1      0.00       N 10.0000
## 2      0.01       N  8.6738
## 3      0.02       N  7.4930
## 4      0.03       N  6.4512
## 5      0.04       N  5.5392
## 6      0.05       N  4.7461
## 9997  49.95       P  0.1635
## 9998  49.96       P  0.1618
## 9999  49.97       P  0.1602
## 10000 49.98       P  0.1586
## 10001 49.99       P  0.1570
## 10002 50.00       P  0.1555

結果：捕食者 (P) と被食者 (N) の個体数動態です。

ggplot(data = dset, aes(x = time, y = value, color=Species)) + geom_line(alpha=0.6,size=2) + labs(x="Time",y="Population size") + theme_bw(20)

plot of chunk unnamed-chunk-16

モデルの状態図 (State plot)

ggplot(data = data.frame(N = dataout$N, P = dataout$P, Model = "Lotka-Volterra 1"), 
    aes(x = N, y = P, color = Model)) + geom_path(alpha = 0.6, size = 2) + geom_point(aes(x = N[1], 
    y = P[1]), size = 4, color = "black") + geom_text(aes(x = N[1], y = P[1]), 
    label = paste("t =", dataout$time[1]), hjust = -0.1, vjust = 0, size = 8, 
    color = "black") + geom_point(aes(x = N[101], y = P[101]), size = 4, color = "black") + 
    geom_text(aes(x = N[101], y = P[101]), label = paste("t =", dataout$time[101]), 
        hjust = -0.1, vjust = 0, size = 8, color = "black") + geom_point(aes(x = N[1001], 
    y = P[1001]), size = 4, color = "black") + geom_text(aes(x = N[1001], y = P[1001]), 
    label = paste("t =", dataout$time[1001]), hjust = -0.1, vjust = 0, size = 8, 
    color = "black") + theme_bw(20)

plot of chunk unnamed-chunk-17

例２。

State = c(N = 10, P = 20)
Param = c(rP = 1, rN = 2, alphaNP = 0.4, alphaPN = 0.1)
Time = seq(0, 50, by = 0.01)
model = ode(State, Time, LV1, Param, method = "ode45")
dataout = as.data.frame(model)

dset = melt(dataout, id="time", variable.name="Species")
dset$Species = as.factor(dset$Species)

Parameters

\( N_{0}= \) 10; \( P_{0}= \) 20; \( r_P= \) 1; \( r_N= \) 2; \( \alpha_{PN}= \) 0.1; \( \alpha_{NP}= \) 0.4

rbind(head(dset), tail(dset))

##        time Species   value
## 1      0.00       N 10.0000
## 2      0.01       N  9.4199
## 3      0.02       N  8.8779
## 4      0.03       N  8.3713
## 5      0.04       N  7.8978
## 6      0.05       N  7.4550
## 9997  49.95       P  0.4140
## 9998  49.96       P  0.4114
## 9999  49.97       P  0.4088
## 10000 49.98       P  0.4062
## 10001 49.99       P  0.4037
## 10002 50.00       P  0.4012

ggplot(data = dset, aes(x = time, y = value, color=Species)) + geom_line(alpha=0.6,size=2) + labs(x="Time",y="Population size") + theme_bw(20)

plot of chunk unnamed-chunk-19 

ggplot(data = data.frame(N = dataout$N, P = dataout$P, Model = "Lotka-Volterra 2"), 
    aes(x = N, y = P, color = Model)) + geom_path(alpha = 0.6, size = 2) + geom_point(aes(x = N[1], 
    y = P[1]), size = 4, color = "black") + geom_text(aes(x = N[1], y = P[1]), 
    label = paste("t =", dataout$time[1]), hjust = -0.1, vjust = 0, size = 8, 
    color = "black") + geom_point(aes(x = N[101], y = P[101]), size = 4, color = "black") + 
    geom_text(aes(x = N[101], y = P[101]), label = paste("t =", dataout$time[101]), 
        hjust = -0.1, vjust = 0, size = 8, color = "black") + geom_point(aes(x = N[1001], 
    y = P[1001]), size = 4, color = "black") + geom_text(aes(x = N[1001], y = P[1001]), 
    label = paste("t =", dataout$time[1001]), hjust = -0.1, vjust = 0, size = 8, 
    color = "black") + ylim(c(0, 22.5)) + theme_bw(20)

plot of chunk unnamed-chunk-20

Lotka-Volterra Model 2

空間や食料がおよぼす競争

つぎの連立微分方程式はロジスティックモデルから構築した。種１ (N1) と種２ (N2) は増加率 \( r \frac{(K-N)}{K} \) によるものです。\( r \) は相対内的増加率であり、生物の最大増加率でもある。\( K \) は環境収容力、すなわちその環境の定員数です。ここまでは、一般的なロジスティック増加モデルですが、あらたに加えた \( \alpha N2 \) や \( \beta N1 \) は相手種との競争率を表している。両種は同じ資源を奪い合っているので、N2が増えるとN1は \( \alpha N2 \) の分減る仕組みになっている。同じようにN2は \( \beta N1 \) の分減少する。 \[ \frac{dN_1}{dt}=r_1 N_1 \left ( \frac{K_1-N_1-\alpha N_2}{K_1} \right ) \\ \frac{dN_2}{dt}=r_2 N_2 \left ( \frac{K_2-N_2-\beta N_1}{K_2} \right ) \]

LV2 = function(Time, State, Param) {
    with(as.list(c(State, Param)), {
        dN1 = r1 * N1 * ((K1 - N1 - alpha * N2)/K2)
        dN2 = r2 * N2 * ((K2 - N2 - beta * N1)/K2)
        return(list(c(dN1, dN2)))
        # The order of this list must match the order of the equations!
    })
}

例１。

State = c(N1 = 10, N2 = 15)
Param = c(r1 = 2, r2 = 1, K1 = 10, K2 = 10, alpha = 0.5, beta = 0.5)
Time = seq(0, 50, by = 0.01)
model = ode(State, Time, LV2, Param, method = "ode45")
dataout = as.data.frame(model)

dset = melt(dataout, id="time", variable.name="Species")
dset$Species = as.factor(dset$Species)

Parameters

\( N1_{0}= \) 10; \( N2_{0}= \) 15; \( r_1= \) 2; \( K_1= \) 10; \( \alpha= \) 0.5; \( r_2= \) 1; \( K_2= \) 10; \( \beta= \) 0.5

データの先頭と末尾を取り出して表示。

rbind(head(dset), tail(dset))

##        time Species  value
## 1      0.00      N1 10.000
## 2      0.01      N1  9.853
## 3      0.02      N1  9.713
## 4      0.03      N1  9.579
## 5      0.04      N1  9.450
## 6      0.05      N1  9.326
## 9997  49.95      N2  6.667
## 9998  49.96      N2  6.667
## 9999  49.97      N2  6.667
## 10000 49.98      N2  6.667
## 10001 49.99      N2  6.667
## 10002 50.00      N2  6.667

結果：Lotka-Volterraの捕食者・被食者モデルと違って、一定の個体数に収束します。

ggplot(data = dset, aes(x = time, y = value, color = Species)) + geom_line(alpha = 0.6, 
    size = 2) + labs(x = "Time", y = "Population size") + theme_bw(20)

plot of chunk unnamed-chunk-24 

ggplot(data = data.frame(N1 = dataout$N1, N2 = dataout$N2, Model = "Lotka-Volterra 2"), 
    aes(x = N1, y = N2, color = Model)) + geom_path(alpha = 0.6, size = 2) + 
    geom_point(aes(x = N1[1], y = N2[1]), size = 4, color = "black") + geom_text(aes(x = N1[1], 
    y = N2[1]), label = paste("t =", dataout$time[1]), hjust = -0.1, vjust = 0, 
    size = 8, color = "black") + geom_point(aes(x = N1[101], y = N2[101]), size = 4, 
    color = "black") + geom_text(aes(x = N1[101], y = N2[101]), label = paste("t =", 
    dataout$time[101]), hjust = -0.1, vjust = 0, size = 8, color = "black") + 
    geom_point(aes(x = N1[1001], y = N2[1001]), size = 4, color = "black") + 
    geom_text(aes(x = N1[1001], y = N2[1001]), label = paste("t =", dataout$time[1001]), 
        hjust = -0.1, vjust = 0, size = 8, color = "black") + ylim(c(4, 17.5)) + 
    xlim(c(5, 12.5)) + theme_bw(20)

plot of chunk unnamed-chunk-25

例2。

State = c(N1 = 10, N2 = 15)
Param = c(r1 = 2, r2 = 1, K1 = 10, K2 = 10, alpha = 1.5, beta = 0.5)
Time = seq(0, 50, by = 0.01)
model = ode(State, Time, LV2, Param, method = "ode45")
dataout = as.data.frame(model)

dset = melt(dataout, id="time", variable.name="Species")
dset$Species = as.factor(dset$Species)

Parameters

\( N1_{0}= \) 10; \( N2_{0}= \) 15; \( r_1= \) 2; \( K_1= \) 10; \( \alpha= \) 1.5; \( r_2= \) 1; \( K_2= \) 10; \( \beta= \) 0.5

結果：この場合、種N2が勝ち抜き、種N1は絶滅する。

ggplot(data = dset, aes(x = time, y = value, color = Species)) + geom_line(alpha = 0.6, 
    size = 2) + labs(x = "Time", y = "Population size") + theme_bw(20)

plot of chunk unnamed-chunk-27 

ggplot(data = data.frame(N1 = dataout$N1, N2 = dataout$N2, Model = "Lotka-Volterra 2"), 
    aes(x = N1, y = N2, color = Model)) + geom_path(alpha = 0.6, size = 2) + 
    geom_point(aes(x = N1[1], y = N2[1]), size = 4, color = "black") + geom_text(aes(x = N1[1], 
    y = N2[1]), label = paste("t =", dataout$time[1]), hjust = 1.1, vjust = 0, 
    size = 8, color = "black") + geom_point(aes(x = N1[51], y = N2[51]), size = 4, 
    color = "black") + geom_text(aes(x = N1[51], y = N2[51]), label = paste("t =", 
    dataout$time[51]), hjust = -0.1, vjust = 0, size = 8, color = "black") + 
    geom_point(aes(x = N1[1001], y = N2[1001]), size = 4, color = "black") + 
    geom_text(aes(x = N1[1001], y = N2[1001]), label = paste("t =", dataout$time[1001]), 
        hjust = -0.1, vjust = 0, size = 8, color = "black") + ylim(c(9, 17.5)) + 
    theme_bw(20)

plot of chunk unnamed-chunk-28

Rosenzweig-MacArthur Model

被食者 (N) の増加率は \( N \times (r_N - \delta N) \) のやなロジスティックモデル から構築したものです。\( r_N \) は出生率そして \( \delta \) は死亡率です。それに、捕食圧率モデル \( \left (\frac{k N}{K+N}P \right ) \) を加える。被食者が増加すると、捕食圧率も増加するが、被食者の個体数が一定の値を超えると、捕食圧は飽和する (\( k \) は最大捕食圧率、\( K \) は最大捕食圧率の半分率を与える被食者の個体数)。捕食圧率モデルと捕食者の個体数を掛け合わせると捕食した被食者の個体数は求められる。さらに、\( \beta \) (変換効率) と掛け合わせると、摂餌による捕食者の個体数の増加も求められる。 \[ \begin{aligned} \frac{dN}{dt} &= N (r_N - \delta N) - \frac{k N}{K + N} P \\ \frac{dP}{dt} &= \beta \frac{k N}{K + N} P - r_P P \end{aligned} \]

RM = function(Time, State, Param) {
    with(as.list(c(State, Param)), {
        dN = rn * N - delta * N^2 - (k * P * N)/(K + N)
        dP = (beta * k * N * P)/(K + N) - rp * P
        return(list(c(dN, dP)))
        # The order of this list must match the order of the equations!
    })
}

State = c(N = 5, P = 2)
Param = c(rn = 1, rp = 1, delta = 1/3, k = 3, K = 1.5, beta = 1)
Time = seq(0, 200, by = 0.01)
model = ode(State, Time, RM, Param, method = "ode45")
dataout = as.data.frame(model)

dset = melt(dataout, id="time", variable.name="Species")
dset$Species = as.factor(dset$Species)

Parameters

\( N_0= \) 5; \( P_0= \) 2; \( r_N= \) 1; \( \delta= \) 0.3333; \( k= \) 3; \( K= \) 1.5; \( \beta= \) 1; \( r_P= \) 1

ggplot(data = dset, aes(x = time, y = value, color = Species)) + geom_line(alpha = 0.6, 
    size = 2) + labs(x = "Time", y = "Population size") + theme_bw(20)

plot of chunk unnamed-chunk-31 

ggplot(data = data.frame(N = dataout$N, P = dataout$P, Model = "Rosenzweig-MacArthur"), 
    aes(x = N, y = P, color = Model)) + geom_path(alpha = 0.6, size = 2) + geom_point(aes(x = N[1], 
    y = P[1]), size = 4, color = "black") + geom_text(aes(x = N[1], y = P[1]), 
    label = paste("t =", dataout$time[1]), hjust = 1.1, vjust = 0, size = 8, 
    color = "black") + geom_point(aes(x = N[101], y = P[101]), size = 4, color = "black") + 
    geom_text(aes(x = N[101], y = P[101]), label = paste("t =", dataout$time[101]), 
        hjust = -0.1, vjust = 0, size = 8, color = "black") + geom_point(aes(x = N[1001], 
    y = P[1001]), size = 4, color = "black") + geom_text(aes(x = N[1001], y = P[1001]), 
    label = paste("t =", dataout$time[1001]), hjust = -0.1, vjust = 0, size = 8, 
    color = "black") + ylim(c(0, 5)) + xlim(c(0, 6)) + theme_bw(20)

plot of chunk unnamed-chunk-32