Overview

In this assignment, you will complete an R script implementing an individual-based simulation of combined ecological (dispersal, reproduction, and death) and evolutionary (adaptive) dynamics. Each individual possesses three attributes: (i) a position on the x dimension, (ii) a position on the y dimension, and (iii) a haploid genotype (0 or 1). Your script will include five functions, but most are written for you (see the code below).

The first function, initPop, creates an initial population on a (by default) 10 x 10 grid. N individuals are generated. Each has a integer-valued location on the x-y grid and a haploid genotype, denoted 0 or 1. The entire population is stored in a list, P, which includes the size of the habitat (size), the number of individuals (N), vectors for the x and y coordinates and the individual genotypes (gen). You do not need to modify this function.

The second function, disperse, allows individuals to disperse (move) to neighboring patches.You will write this function for Task 1 (see details below).

The third function, death, determines which individuals will die at each time step (e.g. year). The version of the function provided below has individuals die with a probability W. You will modify this in Task 2 (below).

The fourth function, reproduce, allows individuals to reproduce. Individuals clone themselves and no new mutations occur (genotype is just copied from the parent). New individuals will have the same x and y coordinates as their parents. The number of offspring per individual is Poisson distributed with mean lambda. You do not need to modify this function.

The fifth and final function, capacity, imposes a carrying capacity, K, for the entire population (across all patches). This is done by selecting K individuals to survive. You do not need to modify this function.

I have also included code (after the function) to run the model and plot the results.

## initialize the population population
initPop <- function(size = 10, N = 20) {
## size specifies number of patches on x and y axis sample
## initial values from 1 to size/2 and -1 to -size/2
  vals <- c(1:(size/2), -1 * 1:(size/2))
  x <- sample(vals, N, replace = TRUE)
  y <- sample(vals, N, replace = TRUE)
  ## sample genotype randomly, 0 or 1
  genotype <- sample(c(0, 1), N, replace = TRUE)
  ## plot the initial population, color denotes genotype
  par(pty = "s")
  plot(x, y, pch = 19, col = c("orange", "blue")[genotype + 1], xlab = "Dimension 1",
      ylab = "Dimension 2", xlim = c(-size/2, size/2), ylim = c(-size/2,
          size/2))
  ## return list, size, N, and x, y and genotype for each
  ## individual
  pop <- list(size = size, N = length(x), x = x, y = y, gen = genotype)
  return(pop)
}

## dispersal function
disperse <- function(P = pop) {
## should allow moves for x and y axis equal probability of move
## to left or right on x and up or down on y
  return(P)
}
## death function, W denotes survival prob.
death <- function(P = pop, W = 0.8) {
## initialize vector to accumulate dead
  dead <- NULL
  ## loop over individuals
  for (i in 1:P$N) {
  ## survive with probability W
    surv <- rbinom(1, 1, prob = W)
    if (surv == 0) {
    ## dead add individual i to dead list
      dead <- c(dead, i)
    }
  }
  ## remove the dead
  Ndead <- length(dead)
  P$N <- P$N - Ndead
  P$x <- P$x[-dead]
  P$y <- P$y[-dead]
  P$gen <- P$gen[-dead]
  if (P$N == 0) {
    print("Popultion extinct")
    }
  return(P)
}
## reproduction, lambda is mean number of offspring per ind.
reproduce <- function(P = pop, lambda = 1.5) {
## empty vectors for new x and y locations of individuals
  xNew <- NULL
  yNew <- NULL
  ## and for genotypes of new individuals
  genotypeNew <- NULL
  ## loop over current individuals
  for (i in 1:P$N) {
  ## sample number of offspring for ind. i
    Nx <- rpois(1, lambda)
    if (Nx > 0) {
    ## add Nx offspring identical to ind i
      xNew <- c(xNew, rep(P$x[i], Nx))
      yNew <- c(yNew, rep(P$y[i], Nx))
      genotypeNew <- c(genotypeNew, rep(P$gen[i], Nx))
    }
  }
  ## determine number of new individuals
  NNew <- length(xNew)
  if (NNew > 0) {
    ## update P by adding new individuals
    P$x <- c(P$x, xNew)
    P$y <- c(P$y, yNew)
    P$gen <- c(P$gen, genotypeNew)
    P$N = P$N + NNew
  }
  return(P)
}

## enforce carrying capacity
capacity <- function(P = pop, K = 10000) {
## if N > K, sample subset to keep
  if (P$N > K) {
    keep <- sample(1:P$N, K, replace = FALSE)
    P$N <- K
    P$x <- P$x[keep]
    P$y <- P$y[keep]
    P$gen <- P$gen[keep]
  }
  return(P)
}

## Run the model, here runs for 100 time steps launch an x11 window
## to see all plots as they are generated x11()
pp <- initPop(20, 50)

for (k in 1:100) {
  pp <- disperse(pp)
  pp <- death(pp)
  pp <- reproduce(pp, 0.4)
  pp <- capacity(pp)
  Sys.sleep(0)
  ## uncomment to plot each step par(pty='s')
  ## plot(jitter(pp$x),jitter(pp$y),pch=19,col=c('orange','blue')[pp$gen+1],,xlab='Dimension
  ## 1', ylab='Dimension
  ## 2',,xlim=c(-pp$size/2,pp$size/2),ylim=c(-pp$size/2,pp$size/2))
}
## plot final, jitter adds some noise to x and y to make all visible
par(pty = "s")
plot(jitter(pp$x), jitter(pp$y), pch = 19, col = c("orange", "blue")[pp$gen +
  1], , xlab = "Dimension 1", ylab = "Dimension 2", xlim = c(-pp$size/2,
  pp$size/2), ylim = c(-pp$size/2, pp$size/2))

Task 1

Your first task is to write the disperse function. This function takes the population as an argument, and returns the updated population. The function should allow each individual to move one square (integer value) left or right and up or down. Specifically, each individual should move left, right or stay in place on the x-axis with equal probability. The same should be true for the y-axis, and moves along the two axes should be independent. I suggest using the sample function for this. See my example below. Moves off of the grid, that is beyond the bounds defined by size, should not be permitted. Try running the simulation before and after adding the disperse function to observe the effect.

## example, using sample to randomly draw -1, 0, or 1 the first
## argument if the vector of values to sample from size gives the
## number of samples to take and replace states whether or not the
## same value can be sampled more than once, which is not relevant
## for size = 1
A <- sample(c(-1, 0, 1), size = 1, replace = FALSE)
A
## [1] -1

Task 2

Your second task is to modify the death function. Specifically, to cause evolution by natural selection, we want to make the probability of death depend on genotype and habitat. Assume the x dimension of the population environment corresponds to some habitat, with different habitat types for positive versus negative values, and assume genotype-dependent survival. In particular, modify the code so that the survival probability for individuals with genotype 0 in negative x-axis squares is 0.8 versus 0.5 in positive x-axis squares. And assume the opposite for individuals with genotype 1, that is a survival probability of 0.8 in positive x-axis squares and 0.5 in negative x-axis squares. Run the model after making these changes to see the effect.

\[ \] \[ \] \[ \] \[ \] \[ \] \[ \]