1 Packages and data

# load the tidyverse package in your current R session
library(tidyverse)

# load the ggplot2 package in your current R session
library(ggplot2)

2 Introduction

The following scripts provide the basic ingredients to run different simulations for each exercise based on the Moran Model. By default, the script runs 20 simulations across 50 years with 100 individuals. The first for-loop analyzes birth and death events, then the second for-loop automates multiple iterations of the simulations over the years.

After creating a plot in each exercise, write a short (1-2 sentences) interpretation of the results.

3 Exercise 1: Random walks to extinction

Run the script below, then plot and interpret the results. This can be achieved in two steps:

  1. Starting with the data frame from the simulation (freq.sp1.df), use gather() to convert your data from wide to long format.
  2. With ggplot, create a plot using geom_line() with “Year” on the x-axis and “Frequency of species 1” on the y-axis. Each simulation (columns in your simulation matrix) is a grouping variable.
# specify the number of simulations (num.sims), the number of years (num.years), and a matrix for output (freq.1.mat)
num.sims <- 20 
num.years <- 50
freq.1.mat <- matrix(ncol = num.sims, nrow = num.years)

# use a for-loop to run through the number of simulations specified in num.sims
# specify parameters and initial conditions of the community
# t0.sp1 represents the frequency of species 1 at time point 0
for (j in 1:num.sims) {
  J <- 100
  t0.sp1 <- 0.5*J  
  community <- vector(length = J)
  community[1:t0.sp1] <- 1
  community[(t0.sp1+1):J] <- 2
  year <- 2
  freq.1.mat[1,j] <- sum(community==1)/J
  
  # second for-loop to run multiple simulations
  # freq.1 represents the frequency of species 1
  # pr.1 represents the frequency of reproduction
  # the last line represents the birth and death rates, based on probabilities
  for (i in 1:(J*(num.years-1))) {
    freq.1 <- sum(community==1)/J
    pr.1 <- freq.1
    community[ceiling(J*runif(1))] <- sample(c(1,2), 1, prob=c(pr.1,1-pr.1))
    
    # record data in the freq.1.mat matrix 
    if (i %% J == 0) {
      freq.1.mat[year, j] <- sum(community==1)/J
      year <- year + 1
    }
  }
}

# set column names in matrix
colnames(freq.1.mat) <- 1:num.sims

# convert freq.1.mat into data frame 
# add a column called year
freq.sp1.df <- as.data.frame(freq.1.mat) %>%
  add_column(year = 1:num.years)

4 Exercise 2: Selective advantage

In the following script, species 1 has a selective advantage over species 2 of 1.2:1 (see “fit.ratio” in the second for-loop).

Plot and interpret the model outcomes using the following scenario: 40 simulations with 200 individuals and 100 years. You will need to change a few numbers in the script below, then plot the results using the same script that you developed in Exercise 1.

# specify the number of simulations (num.sims), the number of years (num.years), and a matrix for output (freq.1.mat)
num.sims <- 20 
num.years <- 50
freq.1.mat <- matrix(ncol = num.sims, nrow = num.years)

# use a for-loop to run through the number of simulations specified in num.sims
# specify parameters and initial conditions of the community
# t0.sp1 represents the frequency of species 1 at time point 0
for (j in 1:num.sims) {
  J <- 100
  t0.sp1 <- 0.5*J  
  community <- vector(length = J)
  community[1:t0.sp1] <- 1
  community[(t0.sp1+1):J] <- 2
  year <- 2
  freq.1.mat[1,j] <- sum(community==1)/J
  
  # second for-loop to run multiple simulations
  # freq.1 represents the frequency of species 1
  # freq.2 represents the frequency of species 2
  # fit.ratio represents the fitness ratio, meaning that species 1 has a selective advantage of 1.2:1
  # pr.1 represents the frequency of reproduction
  # the last line represents the birth and death rates, based on probabilities
  for (i in 1:(J*(num.years-1))) {
    freq.1 <- sum(community==1)/J 
    freq.2 <- 1 - freq.1
    fit.ratio <- 1.2 
    pr.1 <- fit.ratio*freq.1/(fit.ratio*freq.1 + freq.2) 
    community[ceiling(J*runif(1))] <- sample(c(1,2), 1, prob=c(pr.1,1-pr.1)) 
    
    # record data in the freq.1.mat matrix   
    if (i %% J == 0) {
      freq.1.mat[year, j] <- sum(community==1)/J
      year <- year + 1
    }
  }
}

# set column names in matrix
colnames(freq.1.mat) <- 1:num.sims

# convert freq.1.mat into data frame 
# add a column called year
freq.sp1.df <- as.data.frame(freq.1.mat) %>%
  add_column(year = 1:num.years)

5 Exercise 3: Frequency dependent selection

In the following script, we examine how selective advantage and frequency dependence can impact selection.

Plot and interpret the model outcomes using the following scenarios:

5.1 Part I

50 simulations of 500 individuals and 50 years with species 1 having a selective advantage of 1.4

5.2 Part II

20 simulations of 1000 individuals and 30 years with species 1 having a selective disadvantage of 1/1.4. Set a positive frequency dependence of 0.5.

5.3 Part III

50 simulations of 500 individuals and 100 years with species 1 have a selective advantage of 1.2 and a positive frequency dependence of 0.5. Then run the same simulation, but with a negative frequency dependence of -0.5.

# specify the number of simulations (num.sims), the number of years (num.years), and a matrix for output (freq.1.mat)
num.sims <- 20 
num.years <- 50
freq.1.mat <- matrix(ncol = num.sims, nrow = num.years)

# use a for-loop to run through the number of simulations specified in num.sims
# specify parameters and initial conditions of the community
# t0.sp1 represents the frequency of species 1 at time point 0
# fit.ratio.avg represents the average fitness ratio
# freq.dep represents the frequency dependence
for (j in 1:num.sims) {
  J <- 100
  t0.sp1 <- 0.5*J  
  community <- vector(length = J)
  community[1:t0.sp1] <- 1
  community[(t0.sp1+1):J] <- 2
  year <- 2
  fit.ratio.avg <- 1.0
  freq.dep <- 0
  freq.1.mat[1,j] <- sum(community==1)/J
  
  # second for-loop
  # freq.1 represents the frequency of species 1
  # freq.2 represents the frequency of species 2
  # fit.ratio represents the fitness ratio, defined using average fitness ratio and frequency dependence
  # pr.1 represents the frequency of reproduction
  # the last line represents the birth and death rates, based on probabilities
  for (i in 1:(J*(num.years-1))) {
    freq.1 <- sum(community==1)/J
    freq.2 <- 1 - freq.1
    fit.ratio <- exp(freq.dep*(freq.1-0.5) + log(fit.ratio.avg))
    pr.1 <- fit.ratio*freq.1/(fit.ratio*freq.1 + freq.2)
    community[ceiling(J*runif(1))] <- sample(c(1,2), 1, prob=c(pr.1,1-pr.1))
    
    # record data in the freq.1.mat matrix     
    if (i %% J == 0) {
      freq.1.mat[year, j] <- sum(community==1)/J
      year <- year + 1
    }
  }
}

# set column names in matrix
colnames(freq.1.mat) <- 1:num.sims

# convert freq.1.mat into data frame 
# add a column called year
freq.sp1.df <- as.data.frame(freq.1.mat) %>%
  add_column(year = 1:num.years)