Evolution Models

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

knitr::opts_chunk$set(echo = TRUE)

library(ecoevoapps)
library(kableExtra)
library(tidyr)
library(ggplot2)
ggplot2::theme_set(ecoevoapps::theme_apps())

This app presents various models that explore the dynamics between predators and their prey. There are several versions of predator-prey models presented in the tabs below, but they all share a similar structure. The prey species are considered as a single population which, in the absence of the predator, grows either exponentially or logistically. The predator, in turn, depends solely on the prey for its energy (i.e. it is a specialist predator), and if there were no prey available in the system, the predator population simply declines to zero. The more abundant the prey, the more the predator population can grow. But the more the predator population grows, the lower the prey falls. This sets up the cylic dynamics that we see in many predator-prey models.

Historically, versions of the predator-prey models have also been used to study dynamics between parasites and their hosts. Thus, we follow their convention in labeling the "prey" using the variable $H$, and the "predator" using the variable $P$.

{.tabset}

Lotka-Volterra predator-prey model

The classic Lotka-Volterra predator-prey model captures the dynamics between a prey species that has exponential growth, and a predator that consumes prey with no saturation point ("Type 1 functional response", see below):

[ \begin{align} \frac{dH}{dt} &= rH - aHP\ \ \frac{dP}{dt} &= eaHP - dP \end{align} ]

pars_vars <- c("$H$", 
               "$P$", 
               "$r$", 
               "$a$", 
               "$e$", 
               "$d$")
descriptions <- c("Population size of the prey",
                 "Population size of the predator",
                 "Per capita growth rate of the prey",
                 "Attack rate of the predator",
                 "Conversion efficiency of the predator",
                 "Death rate of the predator")
param_df <- data.frame(pars_vars, descriptions)
kable(x = param_df, format = "html", 
      col.names = c("Parameter/Variable", "Description")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed"),
                position = "center")

Equations for the zero net growth isoclines: [ \begin{align} P^ &= \frac{r}{a}\ \ H^ &= \frac{d}{ea} \end{align} ]

sidebarLayout(
  sidebarPanel(
    ### Ask users for parameter values ----
    ## r, a, T_h, e, d

    sliderInput("r_lv_pred1", label = "Per capita growth rate of Prey", 
                min = .0001, max = 1.0, value = .5),
    sliderInput("a_lv_pred1", label = "Predator attack rate", min = .001, 
                max = 1.0, value = .1),
    sliderInput("e_lv_pred1", label = "Predator conversion efficiency",
                min = .001, max = 1.0, value = 0.2),
    sliderInput("d_lv_pred1", label = "Per capita death rate of Predator",
                min = .0001, max = 1.0, value = .3),

    ### Ask users for initial conditions -----
    #N1, N2

    numericInput("H_lv_pred1", label = "Initial population size of Prey", 
                 min = 1, value = 10),
    numericInput("P_lv_pred1", label = "Initial population size of Predator", 
                 min = 1, value = 10),

    ### Ask users for time to simulate ----
    numericInput("t_lv_pred1", label = "Timesteps", min = 10, value = 100),

    checkboxGroupInput("vectors_lv_pred1", label = "Display vector field?", 
                       choices = c("Yes" = "Yes"), selected = "Yes")

  ),  


  mainPanel(renderPlot(plot_lvpred1()),
            renderPlot(np_lvpred1())
            )

)

# Set the initial population sizes
init_lv_pred1 <- reactive({c(H = input$H_lv_pred1, P = input$P_lv_pred1)})

# Set the parameter values
pars_lv_pred1 <- reactive({
  c(r = input$r_lv_pred1, a = input$a_lv_pred1, 
    e = input$e_lv_pred1, d = input$d_lv_pred1)
})

# Time over which to simulate model dynamics
time_lv_pred1 <- reactive({seq(0, input$t_lv_pred1, by = .1)})

# simulate model dynamics
out_lv_pred1 <- reactive({
  data.frame(run_predprey_model(
    time = time_lv_pred1(),
    init = init_lv_pred1(),
    params = pars_lv_pred1()
  ))
})

# Plots ------
## make abundance thru time plot
plot_lvpred1 <- reactive({
  plot_predprey_time(out_lv_pred1())
})

# plot trajectory of population
np_lvpred1 <- reactive({
  if("Yes" %in% input$vectors_lv_pred1) {
    plot_predprey_portrait(out_lv_pred1(), pars_lv_pred1(), 
                                         vectors_field = TRUE)
  } else {
    plot_predprey_portrait(out_lv_pred1(), pars_lv_pred1(), 
                             vectors_field = F)
  }

})

Logistic Prey

This model adds an additional layer of biological realism to the Lotka-Volterra model. Specifically, in the classic model, the prey is modeled as having exponential growth. In this tab, the prey experiences logistic growth up to a carrying capacity $K$.

[ \begin{align} \frac{dH}{dt} &= rH \biggl(1-\frac{H}{K}\biggr) - aHP \ \ \frac{dP}{dt} &= eaHP - dP \end{align} ]

pars_vars <- c("$H$", 
               "$P$", 
               "$r$", 
               "$K$", 
               "$a$", 
               "$e$", 
               "$d$")
descriptions <- c("Population size of the prey",
                 "Population size of predator",
                 "Per capita growth rate of the prey",
                 "Carrying capacity of the prey",
                 "Attack rate of predator",
                 "Conversion efficiency of predator",
                 "Death rate of predator")
param_df <- data.frame(pars_vars, descriptions)
kable(x = param_df, format = "html", 
      col.names = c("Parameter/Variable", "Description")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed"),
                position = "center")

Equations for the zero net growth isoclines: [ \begin{align} P^ &= \frac{r}{a}\bigg(1 - \frac{H}{K}\bigg)\ \ H^ &= \frac{d}{ea} \ \end{align} ]

sidebarLayout(
  sidebarPanel(
    ### Ask users for parameter values ----
    ## r, a, e, d, K

    sliderInput("r_logprey", label = "Per capita growth rate of Prey",
                min = .0001, max = 1.0, value = .5),
    numericInput("K_logprey", label = "Carrying capacity of Prey", 
                 min = 1, value = 100),
    sliderInput("a_logprey", label = "Predator attack rate",
                min = .001, max = 1.0, value = .1),
    sliderInput("e_logprey", label = "Predator conversion efficiency", 
                min = .001, max = 1.0, value = 0.5),
    sliderInput("d_logprey", label = "Per capita death rate of Predator",
                min = .0001, max = 1.0, value = .3),



    ### Ask users for initial conditions -----
    numericInput("H_logprey", label = "Initial population size of Prey",
                 min = 1, value = 10),
    numericInput("P_logprey", label = "Initial population size of Predator", 
                 min = 1, value = 10),

    ### Ask users for time to simulate ----
    numericInput("t_logprey", label = "Timesteps", min = 1, value = 100),

    checkboxGroupInput("vectors_logprey", label = "Display vector field?", 
                       choices = c("Yes" = "Yes"), selected = "Yes")

  ),  


  mainPanel(renderPlot(plot_logprey()), 
            renderPlot(np_logprey())
            )
)

# Set the initial population sizes
init_logprey <- reactive({
  c(H = input$H_logprey, P = input$P_logprey)
})

# Set the parameter values
pars_logprey <- reactive({
  c(r = input$r_logprey, K = input$K_logprey, 
    a = input$a_logprey, e = input$e_logprey, d = input$d_logprey)
  })

# Time over which to simulate model dynamics
time_logprey <- reactive({seq(0, input$t_logprey, by = .1)})

# Simulate the dynamics
out_logprey <- reactive({
    data.frame(run_predprey_model(
    time = time_logprey(),
    init = init_logprey(),
    params = pars_logprey()
  ))
})

# Plots ------

## make abundance thru time plot
plot_logprey <- reactive({
  plot_predprey_time(out_logprey())
})


np_logprey <- reactive({
  if("Yes" %in% input$vectors_logprey) {
    plot_predprey_portrait(out_logprey(), pars_logprey(), 
                                         vectors_field = TRUE)
  } else {
    plot_predprey_portrait(out_logprey(), pars_logprey(), 
                             vectors_field = FALSE)
  }
})

Type II Functional Response

This model goes back to an exponentially growing prey, but adds an additional layer of biological realism to the predator. Specifically, in the classic L-V model, the predator individuals don't spend any time processing they prey, and they never get satiated. In this model, the predator has a Type II functional response, which means that the predator's efficiency in eating the prey is limited by the time it has to invest in processing the prey.

[ \begin{align} \frac{dH}{dt} &= rH - \frac{aH}{1+aT_hH} P\ \ \frac{dP}{dt} &= eP \frac{aH}{1+aT_hH} - dP \end{align} ]

pars_vars <- c("$H$", 
               "$P$", 
               "$r$", 
               "$a$", 
               "$T_{h}$", 
               "$e$", 
               "$d$")
descriptions <- c("Population size of the prey",
                 "Population size of the predator",
                 "Per capita growth rate of the prey",
                 "Attack rate of the predator",
                 "Handling time of the predator",
                 "Conversion efficiency of the predator",
                 "Death rate of the predator")
param_df <- data.frame(pars_vars, descriptions)
kable(x = param_df, format = "html", 
      col.names = c("Parameter/Variable", "Description")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed"),
                position = "center")

Equations for the zero net growth isoclines: [ \begin{align} P^ &= \frac{r}{a} +rT_hH\ \ H^ &= \frac{d}{ea - adT_h} \ \end{align} ]

sidebarLayout(
  sidebarPanel(
    ### Ask users for parameter values ----
    ## r, a, T_h, e, d
    sliderInput("r_lv_pred2", label = "Per capita growth rate of Prey",
                min = .0001, max = 1.0, value = .5),
    sliderInput("a_lv_pred2", label = "Predator attack rate",
                min = .001, max = 1.0, value = .1),
    sliderInput("T_h_lv_pred2", label = "Predator handling time",
                min = 0.001, max = 1.0, value = 0.2),
    sliderInput("e_lv_pred2", label = "Predator conversion efficiency", 
                min = .001, max = 1.0, value = 0.7),
    sliderInput("d_lv_pred2", label = "Per capita death rate of Predator",
                min = .0001, max = 1.0, value = .3),

    ### Ask users for initial conditions -----
    numericInput("H_lv_pred2", label = "Initial population size of Prey",
                 min = 1, value = 10),
    numericInput("P_lv_pred2", label = "Initial population size of Predator",
                 min = 1, value = 10),

    ### Ask users for time to simulate ----    
    numericInput("t_lv_pred2", label = "Timesteps", min = 1, value = 100),

    checkboxGroupInput("vectors_lv_pred2", label = "Display vector field?", 
                       choices = c("Yes" = "Yes"), selected = "Yes")

  ),  


  mainPanel(renderPlot(plot_lvpred2()),
            renderPlot(np_lvpred2())
  )
)

# Set the initial population sizes
init_lv_pred2 <- reactive({
  c(H = input$H_lv_pred2 , P = input$P_lv_pred2)
})

# Set the parameter values
pars_lv_pred2 <- reactive({
  c(r = input$r_lv_pred2, a = input$a_lv_pred2, 
    T_h = input$T_h_lv_pred2 , e = input$e_lv_pred2, 
    d = input$d_lv_pred2 )})

# Time over which to simulate model dynamics
time_lv_pred2 <- reactive({seq(0,input$t_lv_pred2,by = .1)})

# Simulate the dynamics
out_lv_pred2 <- reactive({
  data.frame(run_predprey_model(
    time = time_lv_pred2(),
    init = init_lv_pred2(),
    params = pars_lv_pred2()
  ))
})


# Plots ------

## make abundance thru time plot
plot_lvpred2 <- reactive({
  plot_predprey_time(out_lv_pred2())
})

np_lvpred2 <- reactive({
  if("Yes" %in% input$vectors_lv_pred2) {
    plot_predprey_portrait(out_lv_pred2(), pars_lv_pred2(), 
                                         vectors_field = TRUE)
  } else {
    plot_predprey_portrait(out_lv_pred2(), pars_lv_pred2(), 
                             vectors_field = FALSE)
  }
})

Rosenzweig-MacArthur Model

Finally, Rosenzweig-MacArthur model combines logistic growth for the prey as well as a Type II functional response for the predator.

[ \begin{align} \frac{dH}{dt} &= rH \biggl(1-\frac{H}{K}\biggr) - \frac{aHP}{1+aT_hH} \ \ \frac{dP}{dt} &= e \frac{aHP}{1+aT_hH} - dP \end{align} ]

pars_vars <- c("$H$", 
               "$P$", 
               "$r$", 
               "$K$", 
               "$a$", 
               "$T_{h}$", 
               "$e$", 
               "$d$")
descriptions <- c("Population size of the prey",
                 "Population size of the predator",
                 "Per capita growth rate of the prey",
                 "Carrying capacity of the prey",
                 "Attack rate of the predator",
                 "Handling time of the predator",
                 "Conversion efficiency of the predator",
                 "Death rate of the predator")
param_df <- data.frame(pars_vars, descriptions)
kable(x = param_df, format = "html", 
      col.names = c("Parameter/Variable", "Description")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed"),
                position = "center")

Equations for the zero net growth isoclines: [ \begin{align} P^ &= \frac{r}{a}\bigg(1 - \frac{H}{K}\bigg)(1 + aT_hH)\ \ H^ &= \frac{d}{ea - adT_h} \ \end{align} ]

sidebarLayout(
  sidebarPanel(
    ### Ask users for parameter values ----
    ## r, a, T_h, e, d

    sliderInput("r_rm", label = "Per capita growth rate of Prey", 
                min = .0001, max = 1.0, value = .2),
    numericInput("K_rm", label = "Prey carrying capacity",
                 min = 1, value = 150),
    sliderInput("a_rm", label = "Predator attack rate",
                min = .001, max = 1.0, value = .02),
    sliderInput("T_h_rm", label = "Predator handling time",
                min = 0.001, max = 1.0, value = 0.3),
    sliderInput("e_rm", label = "Predator conversion efficiency", 
                min = .001, max = 1.0, value = 0.6),
    sliderInput("d_rm", label = "Per capita death rate of Predator",
                min = .0001, max = 1.0, value = .4),



    ### Ask users for initial conditions -----
    #N1, N2
    numericInput("H_rm", label = "Initial population size of Prey",
                 min = 1, value = 30),
    numericInput("P_rm", label = "Initial population size of Predator", 
                 min = 1, value = 25),

    ### Ask users for time to simulate ----
    numericInput("t_rm", label = "Timesteps", min = 1, value = 100),

    checkboxGroupInput("vectors_rm", label = "Display vector field?",
                       choices = c("Yes" = "Yes"), selected = "Yes")

  ),

  mainPanel(renderPlot(plot_rm()),
            renderPlot(np_rm())
            )

)

# Set the initial population sizes
init_rm <- reactive({c(H = input$H_rm , P = input$P_rm)})

# Set the parameter values
pars_rm <- reactive({
  c(r = input$r_rm, K = input$K_rm, 
    a = input$a_rm, T_h = input$T_h_rm, 
    e = input$e_rm, d = input$d_rm )})

# Time over which to simulate model dynamics
time_rm <- reactive({seq(0,input$t_rm,by = .1)})

out_rm <- reactive({
    data.frame(run_predprey_model(
    time = time_rm(),
    init = init_rm(),
    params = pars_rm()
  ))
  })

# Plots ------
## make abundance thru time plot
plot_rm <- reactive({
    plot_predprey_time(out_rm())
})

np_rm <- reactive({
  if("Yes" %in% input$vectors_rm) {
    plot_predprey_portrait(out_rm(), pars_rm(), 
                                         vectors_field = TRUE)
  } else {
    plot_predprey_portrait(out_rm(), pars_rm(), 
                             vectors_field = FALSE)
  }
})

Consumer Functional Response {.tabset}

The consumer's functional response captures the relationship between consumer's intake rate of the resource as a function of the resource density. In other words, this is the rate at which consumers can access energy from their prey, as a function of the prey density.

type1_fr <- function (H_init, H_final, a){
  to_return <- c(a*(H_init:H_final))
  return(to_return)
}

type2_fr <- function (H_init, H_final, a, Th){
  to_return <- c((a*(H_init:H_final)/(1+a*Th*(H_init:H_final))))
  return(to_return)
}

Type I

Linear increase in intake rate with resource density, as shown by the equation:

$$f(H) = aH$$

pars_vars <- c("$H$", 
               "$a$")
descriptions <- c("Population size of the prey",
                 "Attack rate of the predator")
param_df <- data.frame(pars_vars, descriptions)
kable(x = param_df, format = "html", 
      col.names = c("Parameter/Variable", "Description")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed"),
                position = "center")

sidebarLayout(
  sidebarPanel(
    ### Ask users for parameter values ----
    # a
    sliderInput("a_fr_type1", label = "Predator attack rate", min = .001, 
                max = 1.0, value = .02),

    ### Ask users for initial and final conditions -----
    #H_init, H_final

    numericInput("H_init_fr_type1", label = "Initial population size of Prey", 
                 min = 1, value = 30),
    numericInput("H_final_fr_type1", label = "Final population size of Prey", 
                 min = 1, value = 100), #Add warning: H_final>H_init 

  ),  

 mainPanel(renderPlot(plot_fr_type1())),

)


# Generate function output and data frame

out_fr_type1 <- reactive({type1_fr(H_init = input$H_init_fr_type1, 
                                   H_final = input$H_final_fr_type1,
                                   a = input$a_fr_type1)})

type1_fr_df <- reactive({data.frame(H = input$H_init_fr_type1:input$H_final_fr_type1,
                                    f_rate = out_fr_type1())})

#Plot

plot_fr_type1 <- reactive ({ggplot(type1_fr_df()) +
                            geom_line(aes(x = H, y=f_rate), size = 2) +
                            scale_color_brewer(palette = "Set1") +
                            xlab("Number of Prey") + 
                            ylab("Foraging Rate") +
                            ylim(0,50)})

Type II

Decelerating intake rate, where the plateau represents consumer saturation:

$$g(H) = \frac {aH}{1+aT_hH}$$

pars_vars <- c("$H$", 
               "$a$", 
               "$T_h$")
descriptions <- c("Population size of the prey",
                 "Attack rate of the predator",
                 "Handling time of the predator")
param_df <- data.frame(pars_vars, descriptions)
kable(x = param_df, format = "html", 
      col.names = c("Parameter/Variable", "Description")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed"),
                position = "center")

sidebarLayout(
  sidebarPanel(
    ### Ask users for parameter values ----
    # a, Th

    sliderInput("a_fr_type2", label = "Predator attack rate", min = .001, 
                max = 1.0, value = .02),
    sliderInput("T_h_fr_type2", label = "Predator handling time", min = 0.001, 
                max = 1.0, value = 0.3),

    ### Ask users for initial and final conditions -----
    #H_init, H_final

    numericInput("H_init_fr_type2", label = "Initial population size of Prey", 
                 min = 1, value = 30),
    numericInput("H_final_fr_type2", label = "Final population size of Prey", 
                 min = 1, value = 100), #Add warning: H_final>H_init 

  ),  

 mainPanel(renderPlot(plot_fr_type2())),

)


# Generate function output and data frame

out_fr_type2 <- reactive({type2_fr(H_init = input$H_init_fr_type2, 
                                   H_final = input$H_final_fr_type2,
                                   a = input$a_fr_type2,
                                   Th = input$T_h_fr_type2)})

type2_fr_df <- reactive({data.frame(H = input$H_init_fr_type2:input$H_final_fr_type2,
                                    f_rate = out_fr_type2())})

#Plot

plot_fr_type2 <- reactive ({ggplot(type2_fr_df()) +
                            geom_line(aes(x = H, y=f_rate), size = 2) +
                            scale_color_brewer(palette = "Set1") +
                            xlab("Number of Prey") + 
                            ylab("Foraging Rate") +
                            ylim(0,4)})

References

Consumer Resource Dynamics by Murdoch, Briggs and Nisbet - digital copy may be available through a library.
Graphical Representation and Stability Conditions of Predator-Prey Interactions, 1963, by M. L. Rosenzweig and R. H. MacArthur.
Primer on functional response and numerical response
Some characteristics of simple types of predation and parasitism, 1959, by C.S. Holling.
Impact of Food and Predation on the Snowshoe Hare Cycle, 1995, by C.J. Krebs et al.
Dynamics of Predation, 2010, by A.N.P. Stevens

Additional references for Host-Parasite/Parasitoid dynamics:
- Overivew of the Nicholson-Bailey model (Links to Wikipedia page)
- Modelling the biological control of insect pests: a review of host-parasitoid models, 1996, by N.J. Mills and W.M. Getz

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

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

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

gauravsk/ecoevoapps
Simulate Dynamics of Ecology/Evolution Models

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

{.tabset}

Lotka-Volterra predator-prey model

Logistic Prey

Type II Functional Response

Rosenzweig-MacArthur Model

Consumer Functional Response {.tabset}

Type I

Type II

References

R Package Documentation

Browse R Packages

We want your feedback!

gauravsk/ecoevoapps Simulate Dynamics of Ecology/Evolution Models

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

{.tabset}

Lotka-Volterra predator-prey model

Logistic Prey

Type II Functional Response

Rosenzweig-MacArthur Model

Consumer Functional Response {.tabset}

Type I

Type II

References

R Package Documentation

Browse R Packages

We want your feedback!

gauravsk/ecoevoapps
Simulate Dynamics of Ecology/Evolution Models