繁殖スケジュールの
リスク分散
議題：繁殖スケジュールのリスク分散は本当か？
空間的リスク分散や時間的リスク分散は、巌佐（1998、「数理生物学入門」）や三浦（2007, 「農業と雑草の生態学」）で簡単な例と簡単な計算でその有利さが示されている。しかし、繁殖スケジュールのリスク分散は下記の議論から有利だと思われるが、定性的あるいは定量的な有利性が上記では紹介されてなかった。本項の目的は、シミュレーションを通じて、繁殖スケジュールのリスク分散の有利性を示すことにある。まずは、繁殖スケジュールとそのリスク分散の紹介をしよう。
生き物には最適な繁殖のタイミングがあるー繁殖スケジュール
いつどのタイミングで、どれくらいの子をなすか、は生物にとって非常に重要な問題である。このような問題は繁殖スケジュールの問題としてよく研究されている。生き物によっては体が大きくなってから子をなすこともあるし、小さい段階で子をなすこともある。体が大きくするのには時間がかかるから、繁殖を開始する体の大きさと繁殖までにかかる時間は相関する。繁殖が遅いほど体が大きくなって1回に産める子も増えるだろうから適応度は高くなるだろう。それなら、ほとんどの生物が体を大きくしてから子をなす、つまり繁殖期が遅くなりそうだが実際はそうではない。例えば、繁殖に不適な環境や極端に死亡率が高くなる環境が訪れる場合はどうだろう。そんな環境が訪れる前に子をなせれば問題ないが、もし自分の体が繁殖期に入ってないのに不適な環境が訪れたら、子をなすことはできずに死亡してしまうかもしれない。ゆえに、たくさんの子をなすために必ずしも繁殖期を遅くする方が適応的とは限らない。一度に産める数は少なくても、繁殖期を速めて子をなす方が適応的な場合もある。
今回は、極端に死亡率が高くなる環境について、雑草学において重要な概念である攪乱をベースに考えてみる。攪乱disturbanceは植物を破壊する物理的な力として定義される（Grime, 1977, American Naturalist）。特に人間活動による除草などがあげられる。当然、繁殖の前に除草されたら子をなすことはできない。
予測不可能な攪乱が繁殖スケジュールに与える影響
植物にとって、人間がいつ草抜きにやってくるかなんて知る由もない。けれども、予測不可能だからと言って、いつまでも繁殖も開始せず、のうのうと生きていたら、除草されて子をなせないだろう。それゆえ、まず予測不可能な攪乱があるときは繁殖期を速めるほうが有利だと考えられる。除草されて1個体も子をなせなかったら元も子もないからだ。次の攪乱が来るまでに、なるべく早くから子をなそうとするとは至極当然な帰結のように見える。
次に、もし攪乱が予測可能ならどうだろう。必ず決まったタイミングで攪乱されるということである。それなら、攪乱される直前まで生育を続けたのちに、最後は体の資源をすべて子に投資して、自分は枯死してしまう、一回繁殖型が有利なことが多い（確か、多回繁殖型と比較できる、引用できる文献があったはず、見つけたら追記）。このような繁殖スケジュールパターンは、季節変動の中でよく観察される。つまり、生育に不適な季節（例えば、冬）は基本的に決まった時期にやってくるから、生育期の最後に花を咲かせて枯死し、種子で休眠してしまう植物はかなり多い。
だが、攪乱の予測が不可能ならどうだろうか。攪乱が比較的、短い期間で起これば、短期間で生育を終えて一回だけ繁殖するほうがわずかでも子をなせて絶滅を免れるかもしれない。攪乱が長く起こらなければ？ 短く生を全うする植物よりも、長く生育して繁殖したほうが子をたくさんなせるだろう。このときに、一気に繁殖までが長い植物に数の差をつけられてしまうかもしれない。しかし、長く生育するリスクは上記で述べたとおりだ。攪乱はいつ起こるかわからないから、一回繁殖型だと繁殖期までが短い場合と長い場合のどちらにとっても高い絶滅リスクにさらされる瞬間が来る。
そこで繁殖期までが短い場合と長い場合のいいとこどり戦略として、無限繁殖型がある（三浦, 2007, 「農業と雑草の生態学」）。雑草ではよくみられる特徴で、生育初期から繁殖を開始し、枯死するまで繁殖を続けるスケジュールである。具体的には、花序が無限に生長を続けて花を咲かせ続ける（無限花序）などを持つことが多い。無限繁殖型は、攪乱が起こるまで体の成長と繁殖を同時に行い、生育期後期になるほど繁殖でなせる子の数が増える。つまり、早くに攪乱が起こっても子をなせ、攪乱が長くおこらなければその分だけたくさん子をなせる。攪乱がいつ起こるかわからないからこそ、特定の繁殖期を持つのではなく、常に繁殖し続ける。ただし、このタイプは、光合成で獲得した資源を体の生長と繁殖に分配し続けるので、一回繁殖型と比べれば生長が遅い（ことを仮定する）。攪乱の予測性が高ければ、一回繁殖型のほうが有利になるのはそのせいである。以上のような、繁殖期の早期化と無限繁殖（無限繁殖の中に繁殖期の早期化を含んでもいいけど）は時間的リスク分散（種子休眠）、空間的リスク分散（種子の長距離散布）と並び、繁殖スケジュールのリスク分散として紹介される（三浦, 2007, 「農業と雑草の生態学」）。
個体群成長速度をシミュレーションする
シミュレーションを行う前に、一回繁殖型と無限繁殖型について、どのような成長をするかを明確に定義する必要がある。下記のような状況を考えることにする。これは、筆者が勝手に「妥当そうな」ものを考えただけである。他の状況になれば今回の定量的議論は結果が変わるので、注意してほしい。あくまで、今回の議論は、一例を示すに過ぎないということである。
下記の図は、一回繁殖型の種1と無限繁殖型の種2の個体の生長速度を示した模式図である。そして、その生長速度の積分（R1とR2）がバイオマスに比例する個体群成長速度を示す。ここでは、単に現在の個体数が次世代では何倍になるか、である。tsは種1の繁殖期の開始、teは種1と種2の生育期の終わりを示す。
上記のように設定した時、B = R1を仮定すると、teは下記のように計算できる。これは蓄積した資源をすべて繁殖に投資することを表している。
また、上記の設定の下でtDにおいて攪乱が起こることを考える。攪乱が起こった後は、その生長速度の積分（R1,2とR2,2）に補正値f（0<f<1）を掛けて、適応度が減少したことを表すことにする。f = 0にしてもいいが、後のシミュレーションで0になるとちょっとめんどくさいことになるので、ここではf > 0とする。
個体群成長速度は攪乱前と攪乱後で合計し、攪乱が起こったタイミングtDに対してR1(tD) = R1,1 + R1,2、R2(tD) = R2,1 + R2,2で表すと、下記のように算出できる（ts、te、a、b、c、fは定数、tDが変数であることに注意）。
tD < 0は攪乱という観点からは植物に影響を与えないはずだが、tD後は適応度が0になるのではなく、下がるとしているので、ここでは芽生えたその瞬間から劣悪な環境、と解釈してもよいだろう。個体群成長速度を算出できたので、あとはこの条件をもとにシミュレーションする。
今回は、a = 2、b = 0.8、c = 0.617、f = 0.001とする。f = 0.001ということは攪乱後の個体群成長速度が1/1000になるということである（つまり実質0）。ts = 1とすれば、上記の式より、te = 1.86となる。fitness関数をtDに応じて下記のように定義する。
------------------------------------------------------
library(plotn)
f <- 0.001
a <- 2
ts <- 1
te <- (2*(exp(a*ts)-1))/(a*exp(a*ts)) + ts
te
## [1] 1.864665
b <- 0.8
c <- 0.617
fitness1 <- function(tD){
if(tD <= ts){
exp(a*ts)*(te-ts)*(f/2)
} else {
if(tD < te){
(1/2)*exp(a*ts)*(f*(te-ts)-(f-1)*(1-(te-tD)/(te-ts))*(tD-ts))
} else {
(1/2)*exp(a*ts)*(te-ts)
}
}
}
fitness2 <- function(tD){
if(tD <= te){
if(tD < 0) tD <- 0
((1-c)/b)*((1-f)*exp(b*tD)+f*exp(b*te)-1)
} else {
((1-c)/b)*(exp(b*te)-1)
}
}
------------------------------------------------------
この条件の下でのtDに応じた個体群成長速度は下記のように図示できる。
------------------------------------------------------
x <- 1.5134388
fitness1(x)
## [1] 1.128458
fitness2(x)
## [1] 1.128458
d1 <- data.frame(x = seq(0,3,length = 200),
sp1 = sapply(seq(0,3,length = 200), fitness1),
sp2 = sapply(seq(0,3,length = 200), fitness2))
plotn(d1[,1], d1[,2:3], mode = "m", type = "l", legend = T,
xlab = "time of disturbance", ylab = "Growth rate",
pt.cex.leg = 0, leg.sp = 4, lty.leg = 1)#図1の図示
overdraw(abline(v = ts, lty = 2),
abline(v = x, lty = 2),
abline(v = te, lty = 2))
------------------------------------------------------
図1 横軸に攪乱が起こるタイミング、縦軸は個体群成長速度を示す。
この図からわかることとして、種1はtD = 1までほとんど個体群成長速度R = 0で、1 < tD < 2までRが上昇する。種2はtD = 2まで一貫してRが増加し続ける。それぞれの様子がちょうど、一回繁殖型と無限繁殖型の様子を表している。なお、tD = 1.51を境に個体群成長速度が種1と2で逆転し、攪乱が起こるまでの時間が長く、生育期を全うできると種1が種2よりも個体群成長速度が高いことを表す。
さて、攪乱の起こるタイミングはランダムに決まるとする。例えば、平均m、標準偏差sの正規分布に従って決まることにする。
まずは、下記のようにm = 2.8、s = 1として、1000世代分の個体群成長速度および個体数を算出する。m = 2.8ということは、tDの分布の平均が2.8であり、teの生育終了から十分離れている。値が大きくなるので対数をとって調整している。forによる繰り返し計算は重たいので、それを避けるために下記のように計算した。目標は下図に示すように、1～1000世代の各世代の個体群成長速度R(tD,g)のベクトルを1000個並べた行列を用意し、その上三角部分を1に置換した行列にする。この行列を、行方向に掛け算すれば1世代目から1～1000世代あとのトータルの個体群成長速度が算出できる。なお、対数をとることにしたので、実際は掛け算ではなく対数の足し算で算出する。
------------------------------------------------------
n1 <- 100 #種1の初期個体数
n2 <- 100 #種2の初期個体数
g <- 1000 #シミュレートする世代数
s <- 1
m <- 2.8
tD <- rnorm(g, m, s) #攪乱タイミングの生成
#以下、1～g世代分の個体群成長速度の対数を算出（対数をとらないと値が極めて大きくなるため）
r_m1 <- matrix(rep(sapply(tD, fitness1), g), ncol = g, byrow = T) * lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v1 <- apply(log10(r_m1), 1, sum)
r_m2 <- matrix(rep(sapply(tD, fitness2), g), ncol = g, byrow = T) * lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v2 <- apply(log10(r_m2), 1, sum)
------------------------------------------------------
これの計算の結果を下記に示す。
------------------------------------------------------
d <- data.frame(G = 0:g, log_n1 = c(log10(n1), log10(n1)+log_r_v1),
log_n2 = c(log10(n2), log10(n2)+log_r_v2)) #1～g世代の個体数の対数
plotn(d[,1], d[,2:3], mode = "m", type = "l", legend = T, ylab = "log(n)", xlab = "Generation") #図2の図示
------------------------------------------------------
図2 m = 2.8。横軸に世代、縦軸は個体数を示す。
この条件の下では、一貫して種1も種2も増加し続ける。だが、種1のほうが増加スピードが速く、個体数で種1は種2を圧倒する。
では、 m = 2.35として、生育期の終了近くで攪乱が起きやすいとどうだろう。
------------------------------------------------------
m <- 2.35
tD <- rnorm(g, m, s) #攪乱タイミングの生成
#以下、1～g世代分の個体群成長速度の対数を算出（対数をとらないと値が極めて大きくなるため）
r_m1 <- matrix(rep(sapply(tD, fitness1), g), ncol = g, byrow = T) * lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v1 <- apply(log10(r_m1), 1, sum)
r_m2 <- matrix(rep(sapply(tD, fitness2), g), ncol = g, byrow = T) * lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v2 <- apply(log10(r_m2), 1, sum)
d <- data.frame(G = 0:g, log_n1 = c(log10(n1), log10(n1)+log_r_v1),
log_n2 = c(log10(n2), log10(n2)+log_r_v2)) #1～g世代の個体数の対数
plotn(d[,1], d[,2:3], mode = "m", type = "l", legend = T, ylab = "log(n)", xlab = "Generation") #図3の図示
------------------------------------------------------
図3 m = 2.35。横軸に世代、縦軸は個体数を示す。
すると、種1の個体数変動が不安定になり始める。が、かろうじてまだ種1の個体数が優勢である。種2は安定して個体数を増殖し続ける。
では、さらにm = 2.22とすると
------------------------------------------------------
m <- 2.22
tD <- rnorm(g, m, s) #攪乱タイミングの生成
#以下、1～g世代分の個体群成長速度の対数を算出（対数をとらないと値が極めて大きくなるため）
r_m1 <- matrix(rep(sapply(tD, fitness1), g), ncol = g, byrow = T) * lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v1 <- apply(log10(r_m1), 1, sum)
r_m2 <- matrix(rep(sapply(tD, fitness2), g), ncol = g, byrow = T) * lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v2 <- apply(log10(r_m2), 1, sum)
d <- data.frame(G = 0:g, log_n1 = c(log10(n1), log10(n1)+log_r_v1),
log_n2 = c(log10(n2), log10(n2)+log_r_v2)) #1～g世代の個体数の対数
plotn(d[,1], d[,2:3], mode = "m", type = "l", legend = T, ylab = "log(n)", xlab = "Generation") #図4の図示
------------------------------------------------------
図4 m = 2.22。横軸に世代、縦軸は個体数を示す。
すると、種1の個体数変動が非常に不安定になる。個体数は対数をとっているので負の値になるということは、個体数は1を下回っており、事実上、種1は絶滅している。だが、まだ種2は安定して個体数を増加できている。
最後、m = 2.05とすると
------------------------------------------------------
m <- 2.05
tD <- rnorm(g, m, s) #攪乱タイミングの生成
#以下、1～g世代分の個体群成長速度の対数を算出（対数をとらないと値が極めて大きくなるため）
r_m1 <- matrix(rep(sapply(tD, fitness1), g), ncol = g, byrow = T) * lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v1 <- apply(log10(r_m1), 1, sum)
r_m2 <- matrix(rep(sapply(tD, fitness2), g), ncol = g, byrow = T) * lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v2 <- apply(log10(r_m2), 1, sum)
d <- data.frame(G = 0:g, log_n1 = c(log10(n1), log10(n1)+log_r_v1),
log_n2 = c(log10(n2), log10(n2)+log_r_v2)) #1～g世代の個体数の対数
plotn(d[,1], d[,2:3], mode = "m", type = "l", legend = T, ylab = "log(n)", xlab = "Generation") #図5の図示
------------------------------------------------------
図5 m = 2.05。横軸に世代、縦軸は個体数を示す。
ここまでくると、種1が一貫して減少し続ける。そして、種2ももはや個体数は増加も減少もしない。
以上の議論から示唆されるのは、生育期の終了から十分たった時に攪乱のタイミングの平均値が存在すれば、種1の一回繁殖型が有利になる。一方で、生育期の終了直後からより短いタイミングに攪乱のタイミングの平均値があると、種2の無限繁殖型が有利になる。次に気になるのは攪乱タイミングが生じるばらつきの影響である。
そこで下記のように平均m、標準偏差sの正規分布から生成された攪乱タイミングから、各世代の個体群成長速度を算出し、さらにその積から、個体群成長速度の幾何平均を求めるgrowth関数を定義する。つまり、この関数は平均m、標準偏差sの正規分布から攪乱タイミングが生成されるときの1世代あたりの平均的な個体群成長速度を算出する。
まず、s = 1.2の時、つまり、攪乱タイミングが非常にばらついて予測ができないとき、
------------------------------------------------------
growth1 <- function(m, n = 10000, sd = 1) 10^(mean(log10(sapply(rnorm(n, m, sd), fitness1))))
growth2 <- function(m, n = 10000, sd = 1) 10^(mean(log10(sapply(rnorm(n, m, sd), fitness2))))
s <- 1.25
d2 <- data.frame(x = seq(0,3,length = 200),
sp1 = sapply(seq(0,3,length = 200), growth1, sd = s),
sp2 = sapply(seq(0,3,length = 200), growth2, sd = s))
plotn(d1[,1], d1[,2:3], mode = "m", type = "l",
legend = T, lty.leg = c(2,2,1,1), pt.cex.leg = 0, leg.sp = 4,
pt.col.leg = c(.default_col[c(1,2)], "blue", "red"),
leg.lab = c("sp1", "sp2", "sp1","sp2"), lty = 2,
xlab = "time of disturbance", ylab = "Growth rate",
mar = c(3.8, 3.8, 2, 1), main = sprintf("sd = %s", s))#図6の図示
overdraw(points(d2[,1], d2[,2], type = "l", col = "blue"),
points(d2[,1], d2[,3], type = "l", col = "red"))
------------------------------------------------------
図6 s = 1.25。横軸に攪乱が起こるタイミング、縦軸は個体群成長速度を示す。点線は攪乱が起こるタイミングがばらつかないとき、実線はばらついているとき。
このとき、攪乱タイミングがばらつかないときに期待できる個体群増殖速度（点線）と比較すると、攪乱タイミングがばらつくとき（実線）は種1と種2の個体群成長速度が逆転するタイミングがより遅くなる。つまり、攪乱が起こるならほとんどのタイミングで種2（無限繁殖型）のほうが種1（一回繁殖型）より有利だということである。
では、順次、ばらつきを小さくする、つまり、攪乱タイミングの予測がしやすくなっていくことを考える。s = 1, 0.7, 0.2の時は下記の通り。
------------------------------------------------------
s <- 1
d2 <- data.frame(x = seq(0,3,length = 200),
sp1 = sapply(seq(0,3,length = 200), growth1, sd = s),
sp2 = sapply(seq(0,3,length = 200), growth2, sd = s))
plotn(d1[,1], d1[,2:3], mode = "m", type = "l",
legend = T, lty.leg = c(2,2,1,1), pt.cex.leg = 0, leg.sp = 4,
pt.col.leg = c(.default_col[c(1,2)], "blue", "red"),
leg.lab = c("sp1", "sp2", "sp1","sp2"), lty = 2,
xlab = "time of disturbance", ylab = "Growth rate",
mar = c(3.8, 3.8, 2, 1), main = sprintf("sd = %s", s))#図7の図示
overdraw(points(d2[,1], d2[,2], type = "l", col = "blue"),
points(d2[,1], d2[,3], type = "l", col = "red"))
s <- 0.75
d2 <- data.frame(x = seq(0,3,length = 200),
sp1 = sapply(seq(0,3,length = 200), growth1, sd = s),
sp2 = sapply(seq(0,3,length = 200), growth2, sd = s))
plotn(d1[,1], d1[,2:3], mode = "m", type = "l",
legend = T, lty.leg = c(2,2,1,1), pt.cex.leg = 0, leg.sp = 4,
pt.col.leg = c(.default_col[c(1,2)], "blue", "red"),
leg.lab = c("sp1", "sp2", "sp1","sp2"), lty = 2,
xlab = "time of disturbance", ylab = "Growth rate",
mar = c(3.8, 3.8, 2, 1), main = sprintf("sd = %s", s))#図8の図示
overdraw(points(d2[,1], d2[,2], type = "l", col = "blue"),
points(d2[,1], d2[,3], type = "l", col = "red"))
s <- 0.2
d2 <- data.frame(x = seq(0,3,length = 200),
sp1 = sapply(seq(0,3,length = 200), growth1, sd = s),
sp2 = sapply(seq(0,3,length = 200), growth2, sd = s))
plotn(d1[,1], d1[,2:3], mode = "m", type = "l",
legend = T, lty.leg = c(2,2,1,1), pt.cex.leg = 0, leg.sp = 4,
pt.col.leg = c(.default_col[c(1,2)], "blue", "red"),
leg.lab = c("sp1", "sp2", "sp1","sp2"), lty = 2,
xlab = "time of disturbance", ylab = "Growth rate",
mar = c(3.8, 3.8, 2, 1), main = sprintf("sd = %s", s))#図9の図示
overdraw(points(d2[,1], d2[,2], type = "l", col = "blue"),
points(d2[,1], d2[,3], type = "l", col = "red"))
------------------------------------------------------
図7 s = 1。横軸に攪乱が起こるタイミング、縦軸は個体群成長速度を示す。点線は攪乱が起こるタイミングがばらつかないとき、実線はばらついているとき。
図8 s = 0.75。横軸に攪乱が起こるタイミング、縦軸は個体群成長速度を示す。点線は攪乱が起こるタイミングがばらつかないとき、実線はばらついているとき。
図9 s = 0.2。横軸に攪乱が起こるタイミング、縦軸は個体群成長速度を示す。点線は攪乱が起こるタイミングがばらつかないとき、実線はばらついているとき。
攪乱タイミングの予測がしやすくなっていくにつれて、種1と種2の個体群成長速度が逆転するタイミングがより早くなる。つまり、攪乱タイミングの予測性が高まるにつれて、種2（無限繁殖型）が種1（一回繁殖型）よりも有利になるタイミングは短くなるということになる。
最後、別の視点から上記と同様の議論を行おう。それは絶滅率である。個体数が1を下回ったら、それは事実上の絶滅である。そこで、攪乱の起こるタイミングが従う正規分布の平均値mと標準偏差sを変えたとき、種1と種2の1000世代分の個体数を計算する。1000世代の中で1度でも個体数が1を下回ったら、そのシミュレーションではその種は絶滅したと判断する。このシミュレーションを100回行って、そのうち何回、絶滅したかを算出する。まずはs = 1.25で攪乱タイミングの予測性が低いとき、
------------------------------------------------------
n1 <- 100
n2 <- 100
g <- 1000 #世代数
re <- 100 #シミュレーションの繰り返し数
v <- c(0, 1.0, 1.5, 1.95, 2.05, 2.22, 2.35, 2.5, 2.65, 2.8, 3)#攪乱タイミングの平均
ex_sp1 <- NULL
ex_sp2 <- NULL
s <- 1.25
for(j in v){
m <- j
ex_sp1_t <- NULL
ex_sp2_t <- NULL
for(i in 1:re){
tD <- rnorm(g, m, s)
r_m1 <- matrix(rep(sapply(tD, fitness1), g), ncol = g, byrow = T)*lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v1 <- apply(log10(r_m1), 1, sum)
r_m2 <- matrix(rep(sapply(tD, fitness2), g), ncol = g, byrow = T)*lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v2 <- apply(log10(r_m2), 1, sum)
log_n1 <- c(log10(n1), log10(n1)+log_r_v1)
log_n2 <- c(log10(n2), log10(n2)+log_r_v2)
ex_sp1_t <- c(ex_sp1_t, sum(log_n1 < 0) > 0) #種1が絶滅したかの判定
ex_sp2_t <- c(ex_sp2_t, sum(log_n2 < 0) > 0) #種2が絶滅したかの判定
}
ex_sp1 <- cbind(ex_sp1, ex_sp1_t) #種1が絶滅したかの判定を各平均につき100回行った結果
ex_sp2 <- cbind(ex_sp2, ex_sp2_t) #種1が絶滅したかの判定を各平均につき100回行った結果
}
d3 <- data.frame(v, sp1 = apply(ex_sp1, 2, function(x) sum(x)/length(x)),
sp2 = apply(ex_sp2, 2, function(x) sum(x)/length(x)))
plotn(d3[,1], d3[,2:3], mode = "m",
ylab = "Extinction rate", xlab = "time of disturbance",
ylim = c(0,1), type = "l", legend = T,
lty.leg = 1, pt.cex.leg = 0, leg.sp = 4,
main = sprintf("sd = %s", s), mar = c(3.8, 3.8, 2, 1)) #図10の図示
------------------------------------------------------
図10 s = 1.25。横軸に攪乱が起こるタイミングの平均、縦軸は絶滅率。
すると、種2（無限繁殖型）のほうが種1（一回繁殖型）に比べてより早いタイミングから、絶滅率が下がり始める。つまり、種1が絶滅してしまうような攪乱タイミングであっても、種2は1000世代生き残れる可能性があるということである。
また、興味深いのが図6を確認すると、種1と種2の個体群成長速度が逆転するのはtD = 2.5付近なのだが、図10を確認すると種1の個体群成長速度が種2より高くなっても種1はまだ十分に絶滅率が高い。これは個体群成長速度のばらつきに依存すると考えられ、平均的に個体群成長速度が高くてもばらつきが大きければ、絶滅する可能性が高いことを表している。したがって、個体群成長速度が種2のほうが小さくても、まだ種2のほうが優先するかの性が残っているということである。
では、s = 1、0.75と攪乱タイミングの予測性が上がった時を確認する。
------------------------------------------------------
ex_sp1 <- NULL
ex_sp2 <- NULL
s <- 1
for(j in v){
m <- j
ex_sp1_t <- NULL
ex_sp2_t <- NULL
for(i in 1:re){
tD <- rnorm(g, m, s)
r_m1 <- matrix(rep(sapply(tD, fitness1), g), ncol = g, byrow = T)*lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v1 <- apply(log10(r_m1), 1, sum)
r_m2 <- matrix(rep(sapply(tD, fitness2), g), ncol = g, byrow = T)*lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v2 <- apply(log10(r_m2), 1, sum)
log_n1 <- c(log10(n1), log10(n1)+log_r_v1)
log_n2 <- c(log10(n2), log10(n2)+log_r_v2)
ex_sp1_t <- c(ex_sp1_t, sum(log_n1 < 0) > 0) #種1が絶滅したかの判定
ex_sp2_t <- c(ex_sp2_t, sum(log_n2 < 0) > 0) #種2が絶滅したかの判定
}
ex_sp1 <- cbind(ex_sp1, ex_sp1_t) #種1が絶滅したかの判定を各平均につき100回行った結果
ex_sp2 <- cbind(ex_sp2, ex_sp2_t) #種1が絶滅したかの判定を各平均につき100回行った結果
}
d3 <- data.frame(v, sp1 = apply(ex_sp1, 2, function(x) sum(x)/length(x)),
sp2 = apply(ex_sp2, 2, function(x) sum(x)/length(x)))
plotn(d3[,1], d3[,2:3], mode = "m",
ylab = "Extinction rate", xlab = "time of disturbance",
ylim = c(0,1), type = "l", legend = T,
lty.leg = 1, pt.cex.leg = 0, leg.sp = 4,
main = sprintf("sd = %s", s), mar = c(3.8, 3.8, 2, 1)) #図11の図示
ex_sp1 <- NULL
ex_sp2 <- NULL
s <- 0.75
for(j in v){
m <- j
ex_sp1_t <- NULL
ex_sp2_t <- NULL
for(i in 1:re){
tD <- rnorm(g, m, s)
r_m1 <- matrix(rep(sapply(tD, fitness1), g), ncol = g, byrow = T)*lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v1 <- apply(log10(r_m1), 1, sum)
r_m2 <- matrix(rep(sapply(tD, fitness2), g), ncol = g, byrow = T)*lower.tri(diag(g), diag = TRUE) + upper.tri(diag(g), diag = F)
log_r_v2 <- apply(log10(r_m2), 1, sum)
log_n1 <- c(log10(n1), log10(n1)+log_r_v1)
log_n2 <- c(log10(n2), log10(n2)+log_r_v2)
ex_sp1_t <- c(ex_sp1_t, sum(log_n1 < 0) > 0) #種1が絶滅したかの判定
ex_sp2_t <- c(ex_sp2_t, sum(log_n2 < 0) > 0) #種2が絶滅したかの判定
}
ex_sp1 <- cbind(ex_sp1, ex_sp1_t) #種1が絶滅したかの判定を各平均につき100回行った結果
ex_sp2 <- cbind(ex_sp2, ex_sp2_t) #種1が絶滅したかの判定を各平均につき100回行った結果
}
d3 <- data.frame(v, sp1 = apply(ex_sp1, 2, function(x) sum(x)/length(x)),
sp2 = apply(ex_sp2, 2, function(x) sum(x)/length(x)))
plotn(d3[,1], d3[,2:3], mode = "m",
ylab = "Extinction rate", xlab = "time of disturbance",
ylim = c(0,1), type = "l", legend = T,
lty.leg = 1, pt.cex.leg = 0, leg.sp = 4,
main = sprintf("sd = %s", s), mar = c(3.8, 3.8, 2, 1)) #図12の図示
------------------------------------------------------
図11 s = 1。横軸に攪乱が起こるタイミングの平均、縦軸は絶滅率。
図12 s = 0.75。横軸に攪乱が起こるタイミングの平均、縦軸は絶滅率。
全ての状況において、種2（無限繁殖型）のほうが種1（一回繁殖型）に比べてより早いタイミングから、絶滅率が下がり始める。また、個体群成長速度が種1のほうがおわ回っていても、種2のほうが絶滅率が低いことも多い。攪乱のタイミングがばらつく、という状況下では、無限繁殖することがいかに絶滅を避けるうえ有効な方法であることがわかるだろう。
おわりに
今回の議論は、個体数増加に関する関数および関数を支配するパラメータをこちらで適当に決めているので、定量的な議論は間違いなく変わるだろう。しかし、定性的には種1と種2の個体群成長速度の関係性が図1のようになっている限り、ほとんど影響を受けないと推測される。このように、無限繁殖型は最大の個体群成長速度が、一回繁殖型に劣っていたとしても、予測不可能な攪乱が起こる環境下では、無限繁殖型のほうがより大きな個体群成長速度を示し、絶滅率が低くなる可能性が高い、と考えられる。