Birth/death 模型随着时间的推移模拟

Question

我正在尝试绘制一个图表，该图表将使用出生死亡模型显示随时间变化的物种数量。到目前为止，我见过的所有包和模型都不允许我输入一个起始数量的物种，这是有问题的，因为我的灭绝率高于我的物种形成率。我对这个问题的数据是；物种形成率（或 b）为 0.0157 消光率（或 d）为 53.3。起始物种数 =250000。目前我已经尝试使用 rbdtree 和 simbdtree。提前谢谢你

Answer 1

假设您不需要生成的物种的系统发育树，而只想知道随着时间的推移物种的数量，那么您真的不需要为此提供一个包：代码相当简单 (如果你知道要实现什么方程）。

同时假设您需要一个随机模拟模型。

为了简单起见，假设您只想知道物种形成或灭绝事件发生时的物种数量（这比弄清楚它们要容易一些在等间隔的时间）。假设人均mutation/speciation/birth比率是b，人均death/extinction比率是d（假设这些比率是恒定的人均是相当标准的，如果你愿意，你可以做出其他假设）。然后，如果当前存在 s 个物种，则事件的总发生率为 (b+d)*s，下一个事件是出生的概率为 b/(b+d)。时间呈指数分布。

num_events <- 1000000
s <- rep(NA, num_events+1)
t <- rep(NA, num_events+1)
t[1] <- 0
s[1] <- 250000
b <- 0.0157
d <- 53.3
set.seed(101)
for (i in 1:num_events) {
   if (s[i]==0) break ## stop if extinct
   delta_t <- rexp(1,rate=(b+d)*s[i])
   t[i+1] <- t[i] + delta_t
   if (runif(1)<(b/(b+d))) s[i+1] <- s[i]+1 else s[i+1] <- s[i]-1
}
plot(t,s,type="s",log="y")
curve(s[1]*exp((b-d)*x), add=TRUE, lwd=3, col=adjustcolor("red", alpha=0.4))

您可以看到，除非您减少到几百个人，否则几乎不值得为随机模拟操心——当数量很大时，理论指数衰减曲线几乎与种群动态匹配。

代码可以很容易地加快一点，但可解释性会付出一些小代价（但最多不会超过几秒钟）。

Birth/death 模型随着时间的推移模拟

Birth/death model over time simulation

r

simulation

phylogeny