#------------------------------------------- # Supplementary R Script accompanying the # manuscript: # # Controlling Invasive Rodents via Synthetic # and the Role of Polyandry # # Manser et al. #------------------------------------------- # Note that this script is based on the most complex # model which allows for polyandry, differential # survival, and mate choice #------------------------------------------- # 1) House keeping #------------------------------------------- rm(list=ls()) # install.packages("deSolve") # download package to solve ODE (if not done yet) library(deSolve) # load package #------------------------------------------- # 2) Write the model function #------------------------------------------- tSryPolyandry <- function(t, state, parameters){ with(as.list(c(state, parameters)), { # Preliminary calculations N = Wx + Wy + Dx + Dy # calculate total population size N y = Dy / (Wy + Dy) # calculate drive frequency among males ft = y * (1 - alpha) / (1 - y*alpha) # calculate the probability of mating with a drive male (Equation S13) p = ft * s * (1 - psi * (1-ft) * (1-r) / (1+r)) # calculate the probability of drive fertilisation (Equation S14) # System of ODE (Equation S12) dWx = b * Wx * (1-p)/2 - (m1 + m2*N)*Wx dWy = b * Wx * (1-p)/2 - (1+s_m) * (m1 + m2*N)*Wy dDx = b * Wx * p/2 - (1+s_m+s_t) * (m1 + m2*N)*Dx dDy = b * Wx * p/2 - (1+s_m+s_t+s_2SRY) * (m1 + m2*N)*Dy + mu # Output return(list(c(dWx,dWy,dDx,dDy))) }) } #------------------------------------------- # 3) Define Parameter Values / Starting Conditions #------------------------------------------- # 3a) Parameters s <- 0.9 # Levels of drive r <- 1 # Levels of drive male sperm competitiveness psi <- 0 # Levels of polyandry alpha <- 0 # Strength of drive male avoidance (mate choice) b <- 3 # Birth rate m1 <- 0.5 # Baseline death rate K <- 1000 # Carrying capacity m2 <- (0.5*b-m1)/K # Density dependent death rate mu <- 0.02*K # Release rate (relative to carrying capacity) mu.start <- 0 # Initial release (relative to carrying capacity) s_t <- 0 # Cost of carrying the t haplotype s_m <- 0 # Cost of being (phenotypically) male s_2SRY <- 0 # Cost of carrying an extra copy of the tSry # 3b) Starting conditions Wx_0 <- K/2 Wy_0 <- K/2 Dx_0 <- 0 Dy_0 <- mu.start*K # 3c) Length of time to be simulated time <- seq(0, 100, length=1000) #------------------------------------------- # 4) Run the model #------------------------------------------- out <- ode(y = c(Wx=Wx_0,Wy=Wy_0,Dx=Dx_0,Dy=Dy_0), times = time, func = tSryPolyandry, parms = c(s=s, r=r, psi=psi, alpha=alpha, b=b, m1=m1, K=K, m2=m2, mu=mu, s_t=s_t, s_m=s_m, s_2SRY=s_2SRY)) #------------------------------------------- # 5) Analyse / Plot Model Output #------------------------------------------- # 5a) Calculations N <- out[,'Wx']+out[,'Wy']+out[,'Dy']+out[,'Dx'] # Calculate Population Sizes y <- out[,'Dy'] / (out[,'Dy']+out[,'Wy']) # Calculate Drive Frequency (among fertile males) time <- out[,'time'] # 5b) Plot opacity <- 0.8 plot(time,N,ylim=c(0,max(N,na.rm=T)),type="n",ylab="Density",xlab="Time",axes=F) polygon(c(time,rev(time)), c(out[,'Wx'],rep(0,length=length(time))) ,col=rgb(1,0,0,opacity),border=NA) polygon(c(time,rev(time)), c(out[,'Wx']+out[,'Wy'],rev(out[,'Wx'])) ,col=rgb(0,1,0,opacity),border=NA) polygon(c(time,rev(time)), c(out[,'Wx']+out[,'Wy'],rev(out[,'Wx']+out[,'Wy']+out[,'Dy'])),col=rgb(0,0,1,opacity),border=NA) polygon(c(time,rev(time)), c(N,rev(out[,'Wx']+out[,'Wy']+out[,'Dy'])) ,col=rgb(1,1,0,opacity),border=NA) axis(1);axis(2)