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

knitr::opts_chunk$set(echo = TRUE)
library(tidyr)
library(ggplot2)
library(purrr)
library(deSolve)
library(ecoevoapps)
library(patchwork)
library(knitr)
library(kableExtra)

dndt <- function(params) {
  r <- params["r"]
  if(!is.na(params["K"])) {
    K <- params["K"]
    (1:K)*r*(1-(1:K)/K)
  } else {
    r*1:10000
  }
}

dnNdt <- function(params) {
  r <- params["r"]
  if(!is.na(params["K"])) {
    K <- params["K"]
    r*(1-seq(1:K)/K)
  } else {
    rep(r,10000)
  }
}

This is an entry in the EcoEvoApps project, a collection of freely available apps to simulate canonical models from theoretical ecology and evolution. For more apps in this collection, please visit the project website: https://ecoevoapps.gitlab.io.

Model description and interactive app

{.tabset}

English

Turkish

中文

Model dynamics

{.tabset}

Interactive application

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.01, max = 1, value = .2, step = 0.01),

    # 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.01, max = 1, value = .15, step = 0.01),

                     # 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), size = 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), size = 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), size = 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), size = 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, max(input$r_val, input$r_val2)*1.2), 
                         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')
})

Simulating model dynamics in R

Further reading

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

"Exponential & logistic growth" at Khan Academy.

---

App maintained by Gaurav Kandlikar ([e-mail])

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