In gauravsk/ecoevoapps: Simulate Dynamics of Ecology/Evolution Models

中文 | Español | English | português | Turkish

source("single_populations_app_setup.R")

本文介绍单个物种的种群动态。为了方便起见，我们只考虑无个体迁入或迁出的封闭种群。该种群大小因个体出生而增长，因个体死亡而缩小。若种群增长不受抑制，则其增长呈指数型。无论种群大小变化，其出生率和死亡率均不变。已知人均出生率（per-capita birth rate）为$b$， 人均死亡率（per-capita death rate）为$d$， 则种群大小$N$的变化为：

$$\frac{dN}{dT} = bN - dN = (b-d)N$$

$(b-d)$ 是一个常数，可合并写为自然增长率(intrinsic growth rate)$r$。则指数型增长可写为：

$$\frac{dN}{dT} = rN$$

对于一个指数型增长的种群，其种群大小随时间变化的轨迹呈J型曲线。而J的形状由$r$的数值决定。👇请使用文末的app探索$r$值如何改变种群的轨迹。

限制种群增长

上文的指数增长模型（exponential growth model）预测，该类种群的长期规模会趋于无穷大。但现实中，没有种群的数量能够无限制增长，因此指数增长模型可以被改进。

指数增长模型的前提条件是人均出生率$b$和人均死亡率$d$均不变。我们可以放宽这个条件，并假设当种群的增长趋于环境承载力(carrying capacity $K$)时，其总增长率(net growth rate $rD$)呈线性递减。此类种群的动态变化可用逻辑斯谛增长(logistic growth)模型描述：

$$\frac{dN}{dt} = rN \left(1-\frac{N}{K}\right)$$

呈逻辑斯谛增长的种群的会持续增长，直到种群规模与环境承载力相等($N = K$)。当种群规模超出环境承载力时，种群出现负增长，直到其规模减少至环境承载力。

滞后逻辑斯谛增长

在某些情况下，一个种群的增长率会取决于该种群在过去某个时间点的规模。比如，生物的繁殖产量会与其未成年时资源竞争的强度相关。对于这种情况，某时间点$t$时的种群增长率可以用过去$t-\tau$时的种群大小表示：

$$ \frac{dN}{dt} = rN_{t} \left(1-\frac{N_{t-\tau}}{K}\right)$$ 在app中可选择密度制约，并请勾选“Lagged density dependence”（滞后的密度依赖）选项，来模拟滞后逻辑斯谛增长。

参数表

pars_vars <- c("$r_i$", 
               "$K_i$", 
               "$N_i$",
               "$\\tau$")
descriptions <- c("物种$i$的自然增长率",
                 "物种$i$的环境承载力（carrying capacity）",
                 "物种$i$的种群大小",
                 "滞后逻辑斯谛增长的滞后时间")
param_df <- data.frame(pars_vars, descriptions)
kable(x = param_df, format = "html", 
      col.names = c("参数/变量", "描述")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed"),
                position = "center")

Interactive app

sidebarLayout(
  sidebarPanel(
    # User defined density dependent or independent growth ----
    radioButtons("densdep", "Type of population growth", 
                 choices = c(`Density Independent` = "indep", `Density Dependent` = "dep")),

    # User defined intrinsic growth rate r -----
    sliderInput("r_val", label = "Intrinsic growth rate (r):",
                min = -0.25, max = 1, value = .25),

    # If Density Dependent growth requested, ask for K -----
    conditionalPanel(condition = "input.densdep == 'dep'",
                     sliderInput("K_val", label = "Carrying capacity (K)",
                                 min = 1, max = 10000, value = 500, step = 10)),
    # If Density Dependent growth requested, ask if time lag requested -----
    conditionalPanel(condition = "input.densdep == 'dep'",
                     radioButtons("densdepLag", label = "Delayed density dependence?",
                                  choices = c(`No lag` = "nolag",`Lagged density dependence` = "lag"))),

    # If timelag requested in density-dependent growth, ask for tau -----
    conditionalPanel(condition = "input.densdepLag == 'lag'",
                     sliderInput("tau", label = "Time lag (tau)",
                                 min = 0, max = 4, value = 1, step = .1)),

    # User defined initial population size  -----
    sliderInput("N1_val", label = "Initial population size:",
                min = 1, max = 1000, value = 50, step = 10),
    hr(),

    # User can add a second species ----
    checkboxInput("secondsp", "Add a second species?", value = F), 

    # If second species requested, add the following... ----
    conditionalPanel(condition = "input.secondsp == true",
                     # User defined density dependent or independent growth for sp2 ----
                     radioButtons("densdep2", "Type of population growth for Species 2", 
                                  choices = c(`Density Independent` = "indep", `Density Dependent` = "dep")),

                     # User defined intrinsic growth rate r2 -----
                     sliderInput("r_val2", label = "Intrinsic growth rate for species 2 (r2):",
                min = -0.25, max = 1, value = .25),

                     # If Density Dependent growth requested for sp2, ask for K2 -----
                     conditionalPanel(condition = "input.densdep2 == 'dep'",
                                      sliderInput("K_val2", label = "Carrying capacity for species 2 (K2)",
                                                  min = 1, max = 10000, value = 500, step = 10),
                     ),
                     # If Density Dependent growth requested, ask if time lag requested -----
                     conditionalPanel(condition = "input.densdep2 == 'dep'",
                                      radioButtons("densdepLag2", label = "Delayed density dependence in Species 2?",
                                                   choices = c(`No lag` = "nolag",`Lagged density dependence` = "lag"))),

                     # If timelag requested in density-dependent growth, ask for tau -----
                     conditionalPanel(condition = "input.densdepLag2 == 'lag'",
                                      sliderInput("tau2", label = "Time lag sp 2 (tau2)",
                                                  min = 0, max = 4, value = 1, step = .1)),
                     sliderInput("N2_val", label = "Initial population size of species 2:",
                                 min = 1, max = 1000, value = 50, step = 10),
                         hr()),

    # User can select which plots to show ----
    checkboxGroupInput("plots_to_show", "Select which plots to show",
                       c("N vs. Time" = "Nvtime",
                         "ln(N) vs. Time" = "lnNvtime",
                         "dN/dt vs. N" = "n1overn",
                         "dN/Ndt vs. N" = "n1novern"), selected = c("Nvtime", "lnNvtime")),
    hr(),

    # User defined length of simulation -----
    numericInput("max_time", label = "Length of simulation:",
                value = 100),

                     # User defined initial population size  -----
    ),



  # Panel of plots -----
  mainPanel(
    renderPlot({plots_to_print()}, width = 600,
               height = 600)
  )
)

# Get user defined parameters for sp 1  ------
params <- reactive({
  if(input$densdep == "indep") {
    c(r = input$r_val) 
  } else if(input$densdepLag == 'lag') {
    c(r = input$r_val, K = input$K_val, tau = input$tau) 
  } else {
    c(r = input$r_val, K = input$K_val) 
  }
})
init <- reactive({c(N1 = input$N1_val)})
time <- reactive({seq(from = 0, to = input$max_time, by=1)})

params2 <- reactive({
  if(input$densdep2 == "indep") {
    c(r = input$r_val2) 
  } else if(input$densdepLag2 == 'lag') {
    c(r = input$r_val2, K = input$K_val2, tau = input$tau2) 
  } else {
    c(r = input$r_val2, K = input$K_val2) 
  }
})
init2 <- reactive({c(N1 = input$N2_val)})

# Generate trajectories for sp 1 --------
pop_growth <- reactive({
  if(input$densdep == 'indep') {
    data.frame(ecoevoapps::run_exponential_model(time= time(), init = init(), params = params()))
  } else {
    data.frame(ecoevoapps::run_logistic_model(time= time(), init = init(), params = params()))
  }
})

# Generate trajectories for sp 2 --------
pop_growth2 <- reactive({
  if(input$secondsp) {
    if(input$densdep2 == 'indep') {
      data.frame(ecoevoapps::run_exponential_model(time= time(), init = init2(), params = params2()))
    } else {
      data.frame(ecoevoapps::run_logistic_model(time= time(), init = init2(), params = params2()))
    }
  }
}) 

pops_combined <- reactive({
  if(input$secondsp) {
    dplyr::bind_rows(pop_growth(), pop_growth2(),
                     .id = "which_pop")
  } else {
    pg <- pop_growth()
    pg$which_pop = "1"
    pg
  }
})
# Make plots -----------

# 0. Make caption  ------
plot_caption1 <- reactive({
  if(input$densdep == 'indep') {
  paste0("Species 1 parameter values: r1 = ", input$r_val)
  } else if (input$densdepLag == "nolag") { 
    paste0("Species 1 parameter values: r1 = ", input$r_val, ", K1 = ", input$K_val) 
    } else { 
    paste0("Species 1 parameter values: r1 = ", input$r_val, ", K1 = ", input$K_val, ", Tau1 = ", input$tau) 
    }
})

plot_caption2 <- reactive({
  if(input$secondsp) {
    if(input$densdep2 == 'indep') {
      paste0("Species 2 parameter values: r2 = ", input$r_val2)
    } else if (input$densdepLag2 == "nolag") { 
      paste0("Species 2 parameter values: r2 = ", input$r_val2, ", K2 = ", input$K_val2) 
    } else { 
      paste0("Species 2 parameter values: r2 = ", input$r_val2, ", K2 = ", input$K_val2, ", Tau2 = ", input$tau2) 
    }
  }
})

plot_caption <- reactive ({
  if(input$secondsp) {
    paste0(plot_caption1(), "\n", plot_caption2())
  } else {
    plot_caption1()
  }
})

# 1. Trajectory of N vs. time -----------
trajectory <- reactive({
  if("Nvtime" %in% input$plots_to_show) {
    ggplot(pops_combined()) + 
      geom_line(aes(x = time, y = N1, color = which_pop), linewidth = 2) +
      ecoevoapps::theme_apps() +
      scale_x_continuous(expand = c(0, 0, .1, 0)) +
      scale_y_continuous(expand = c(0, 0, .1, 0)) +
      annotate("text", x = 0, y = 0, label = "") +
      scale_color_brewer("Species", palette = "Set1") +
      ylab("Population size") 
    } 
})

# 2. Trjectory of ln(N) vs. time --------
trajectory_log <- reactive({
  if("lnNvtime" %in% input$plots_to_show) {
    ggplot(pops_combined()) + 
      geom_line(aes(x = time, y = log(N1), color = which_pop), linewidth = 2) +
      scale_color_brewer("Species", palette = "Set1") +
      ecoevoapps::theme_apps() +
      scale_x_continuous(expand = c(0, 0, .1, 0)) +
      scale_y_continuous(expand = c(0, 0, .1, 0)) +
      annotate("text", x = 0, y = 0, label = "") +
      ylab("ln(Population size)") 
    }
})

# 3. dN/dT vs. N -----------
n1_over_n <- reactive({
  if("n1overn" %in% input$plots_to_show) {

    dndt_out <- data.frame(dndt = dndt(params()), which_pop = "1")
    dndt_out$x <- 1:nrow(dndt_out)

    dndt_out_for_plot <-
      if(input$secondsp) {
        dn2dt_out <- data.frame(dndt = dndt(params2()), which_pop = "2")
        dn2dt_out$x <- 1:nrow(dn2dt_out)
        dplyr::bind_rows(dndt_out, dn2dt_out)
      } else {
        dndt_out
      }

    ggplot(dndt_out_for_plot) +
      geom_line(aes(x= x, y = dndt, color = which_pop), linewidth = 2) + 
      xlab("Population Size") + 
      ylab("Population growth rate (dN/dt)") + 
      scale_color_brewer("Species", palette = "Set1") +
      scale_x_continuous(expand = c(0, 0, .1, 0)) +
      scale_y_continuous(expand = c(0, 0, .1, 0)) +
      theme_apps() 
    }
})

# 4. dN/Ndt vs. N -----------
n1n_over_n <- reactive({
  if("n1novern" %in% input$plots_to_show) {

    dnNdt_out <- data.frame(dnNdt = dnNdt(params()), which_pop = "1")
    dnNdt_out$x <- 1:nrow(dnNdt_out)

    dnNdt_out_for_plot <-
      if(input$secondsp) {
        dnN2dt_out <- data.frame(dnNdt = dnNdt(params2()), which_pop = "2")
        dnN2dt_out$x <- 1:nrow(dnN2dt_out)
        dplyr::bind_rows(dnNdt_out, dnN2dt_out)
      } else {
        dnNdt_out
      }

    ggplot(dnNdt_out_for_plot) +
      geom_line(aes(x = x, y = dnNdt, color = which_pop), linewidth = 2) + 
      xlab("Population Size") + 
      ylab("per-capita population growth rate\n(dN/Ndt)") + 
      scale_color_brewer("Species", palette = "Set1") +
      scale_x_continuous(expand = c(0, 0, .1, 0)) +
      # limits = c(0, min(max(pop_growth()$N1)), max(pop_growth2()$N1)
      scale_y_continuous(limits = c(-0.26, 1),                           expand = c(0, 0, .1, 0)) +
      theme_apps() 
    }
})

# Make a list of plots and print out plots based on which ones were requested ----
plot_list <- reactive({
  list(trajectory(), trajectory_log(),
       n1_over_n(), n1n_over_n()) %>% 
    discard(is.null)
  })

plots_to_print <- reactive({
  wrap_plots(plot_list(), ncol = 2) + 
    labs(caption = plot_caption()) + 
    plot_layout(guides = "collect") & theme(legend.position = 'bottom')
})

Further reading

"How populations grow: the Exponential and the logistic", John Vandermeer, Nature Education Knowledge.

"Exponential & logistic growth" at Khan Academy.

---

suppressWarnings(ecoevoapps::print_app_footer(language = "ch"))

gauravsk/ecoevoapps documentation built on July 9, 2024, 9:37 p.m.