#------------------------------------------------------------------------------ ## MetaPop #------------------------------------------------------------------------------ # by Mandy Barron # version 6.1, 26 April 2018 # edited by John Kean # version 7, 9 May 2019 # Simulates early establishment and spread of an invasive species # in a landscape of habitat patches. # Simulates effects of biopesticide applications and/or habitat removal. #------------------------------------------------------------------------------ # Preliminaries # clean up variables rm(list=ls()) gc() #set working directory to source file location dirN<-dirname(rstudioapi::getSourceEditorContext()$path) setwd(dirN) #------------------------------------------------------------------------------ # Functions # randomPatch # generates habitat map with desired % habitat and level of fragmentation # simmap.R returns habitat map as a vector (c.f. simmap2.R which returns as raster) randomPatch <- dget(paste(dirN,"/simmap.R",sep='')) # roundUp # rounds input up to the nearest multiple of "to" roundUp <- function(x, to=10) { to*(x%/%to + as.logical(x%%to)) } # genRegPatchLocs # calcs xy's for each habitat patch as the cell centres genRegPatchLocs <- function(x,y,griddim) { st<-round(griddim/2,1) endx<-x-round(griddim/2,1) endy<-y-round(griddim/2,1) xys <- expand.grid(x=seq(st,endx,by=griddim), y=seq(st,endy, by=griddim))+0 return(xys) } # numSurv # returns integer subset surviving according to a Ricker function # Nin = input population size # Kin = switch: if Kin = 0 then nothing survives # diSurv = density independent survival proportion # dd = density dependent survival rate numSurv <- function(Nin,Kin,diSurv,dd=0) { surv <- ifelse(Kin>0, diSurv*exp(-dd*Nin), 0) ifelse(Nin<10000, rbinom(length(Nin), size=Nin, prob=surv), round(Nin*surv)) } # genWrapdDistMat # Calculates distance between grid cells centres but on a torus so edges # considered adjacent i.e. cell(2,1) is considered next to cell (2,60) genWrapdDistMat<-function(gridCoords,lx,ly) { halfwidth<-lx/2 halfht<-ly/2 distMat <- matrix(0,nrow=nrow(gridCoords),ncol=nrow(gridCoords)) for(ci in 1:nrow(gridCoords)) for (ri in 1:nrow(gridCoords)) { dx<-abs(gridCoords$x[ci]-gridCoords$x[ri]) dx<-ifelse(dx>halfwidth,lx-dx,dx) dy<-abs(gridCoords$y[ci]-gridCoords$y[ri]) dy<-ifelse(dy>halfwidth,ly-dy,dy) distMat[ri,ci]<-sqrt((dx^2)+(dy^2)) } return(distMat) } # probDisp # calculates dispersal probability according to half-normal distribution # and sets probs to zero for those outside maximum distance so don't deal with tiny numbers probDisp<-function(Dists,dispScal,maxDist) { pd<-(1/lpdnorm)*exp(-(Dists^2)/(2*dispScal^2)) pd[which(Dists>maxDist)]<-0 return(pd) } # probMated # calculates probability of female being mated (once or more) from number of males in # and the distance to surrounding cells probMated<-function(Males,Dists,distScal,maxProb) { notmated<-prod((1-maxProb*exp(-distScal*Dists))^Males) 1-notmated } # numProd # calculates mean number of eggs laid per cell according to number of mated females # (can be density-dependent although we didn't use this) then draws actual number # from a Poisson dist.. numProd<-function(Nffem,Kin,diProd,dd=0) { produc<-ifelse(Kin>0,Nffem*diProd*exp(-dd*Nffem/Kin),0) ifelse (produc<1000,rpois(length(Nffem),lambda = produc),produc) } # sseq # calculates which cells to treat for different swath widths sseq <- function(from, to, by){ times <- (to-from) %/% sum(by) x <- cumsum(c(from, rep(by, times+1))) x[x<=to] } #------------------------------------------------------------------------------ # Parameters # Landscape parameters landscape.x<-1200 # metres wide landscape.y<-1200 # metres long gridcell<-20 # patch size in metres along each square side # Population growth parameters meanK<-250000 # mean carrying capacity per ha corresponding to 'outbreak' levels of fourth instar larvae; Campbell, 1981 sexRatio<-0.5 # proportion eggs laid which are female eggsPerFem<-400 # eggs laid per female EtoLsurv<-0.5 # proportional survival of eggs LtoAsurv<-0.1 # maximum proportional survival of larvae and pupae # Dispersal parameters Ldiffcoeff<-0.003 # larval diffusion coefficient in km2 from Robinet et al. 2008 MattrDistcoef<-0.01 # distance scalar for males attracted to females, 0.01 is the mean of those in Robinet et al. 2008 MattrMult<-0.01 # maximum probability of male locating a female estimated by simulation to match Figure 2 of Robinet at al. 2008 # Simulation parameters and initial conditions tMax<-10 # length of each simulation in years numEggMasses<-5 # initial number of egg masses initNumPatchOcc<-1 # initial munber of patches occupied numReps<-500 # number of replicate simulations # Calculated parameters # Landscape parameters nx<-landscape.x/gridcell # number of cells along x direction ny<-landscape.y/gridcell # number of cells along y direction cell.area<-gridcell^2 # area of each habitat patch in m2 landscape.area<-landscape.x*landscape.y # area of simulation arena in m2 numPatch<-(landscape.x*landscape.y)/gridcell^2 # number of habitat patches # Population parameters meanKperCell<-meanK*cell.area/10000 # carrying capacity per patch rmax<-log(eggsPerFem*EtoLsurv*LtoAsurv) # maximum rate of increase LtoAddcoef<-rmax/meanKperCell # density dependence coefficient for larval-pupal stage transition # Dispersal parameters LdispDistScal<-round(sqrt(2*Ldiffcoeff*1)*1000,2) #in metres, derived as: sqrt(2*D*t) where D is the diffusion coefficient # calc maximum distance to evaluate dispersal over (5*SD) and calculate the # normalising constant (takes a lot of processing to evaluate this numerically!) maxLdispRad<-5*LdispDistScal #maximum larval dispersal radius in m ncells<-roundUp(maxLdispRad,gridcell)/gridcell #maximum larval dispersal radius in cells xy <- expand.grid(seq(-(ncells*gridcell),ncells*gridcell,by=gridcell), seq(-(ncells*gridcell),ncells*gridcell,by=gridcell))+0 dmat<-(0-xy[,1])^2+(0-xy[,2])^2 lpdmat<-exp(-dmat/(2*LdispDistScal^2)) lpdnorm<-sum(lpdmat) rm(list = c("lpdmat","dmat","xy")) gc() #------------------------------------------------------------------------------ # Management interventions # Sets up a "control schedule" ctrlTypes<-c("Pesti","NatEnms","PherTrap","HabDestr") ctrlStages<-c("E","L","M","Fe") ctrlPatches<-c(1:numPatch) #create a dummy control schedule (stage, type, intensity, what patches to control, when to control) control<-list(stage=ctrlStages[1],type=ctrlTypes[1],intensity=0.9,patches=ctrlPatches[1:3],times=rep(0,times=tMax)) #------------------------------------------------------------------------------ # Lists of parameter combinations to simulate # comment in or out for different sets of simulations # swath width (in cells) for habitat removal or pesticide application swathWlist<-c(1,2,3,5,6,10,15,30) # proportion of habitat to randomly remove #pHabRemlist<-seq(from=0.70, to=0.86, by=0.02) # proportional pesticide efficacy (% kill) #pPestEfflist<-seq(from=0.80, to=0.96, by=0.02) #initEMlist<-1:10 #eggsPEMlist<-c(200,300,400,500) #propHablist<-seq(from=0.1, to=1, by=0.1) #fragParmlist<-seq(from=0.05, to=5, by=0.05) #------------------------------------------------------------------------------ # Initialise the grid # can take all this outside tmts and reps loop if all reps have the same landscape patchxys<-genRegPatchLocs(landscape.x,landscape.y,gridcell) # use these distances between patches for absorbing boundaries #distBWpatches<-as.matrix(dist(patchxys, method = "euclidean", diag = TRUE, upper=TRUE)) # use these distances between patches for wrapped/revolving boundaries distBWpatches<-genWrapdDistMat(patchxys,landscape.x,landscape.y) #------------------------------------------------------------------------------ # Outer loop for simulations of different parameter combinations for (sw in 1:length(swathWlist)) { swathwidth<-swathWlist[sw] # create a multi-dimensional array to hold results for this parameter combination # dimensions = time*3statevars*numreps metaResArr<-array(data = NA, dim = c((tMax+1),3,numReps), dimnames = list(Time=seq(0,tMax,1),State=c("sumFe","sumEggs","numOccEggs"),repl=seq(1,numReps,1))) # ceate a control schedule # if want to randomly allocate patches to control need to do this inside the replicate loop # otherwise for fixed swath widths can allocate patches to control ($patches) here once for all replicates # here ctrl1 is a one off habitat destruction implemented by setting cell carrying capacity to zero ctrl1<-control ctrl1$stage<-ctrlStages[1] ctrl1$type<-ctrlTypes[4] ctrl1$intensity<-1 ctrl1$times[4]<-1 # tartan pattern same width removed as left-> 0.5*0.5 = 25% remaining as habitat = 75% cells removed swincrement<-(swathwidth+1)*gridcell st<-round(gridcell/2,1) endx<-landscape.x-round(gridcell/2,1) endy<-landscape.y-round(gridcell/2,1) byset<-c(rep(gridcell,times=swathwidth-1),swincrement) xset<-sseq(st, endx, by=byset) yset<-sseq(st, endy, by=byset) patchIDs<-c(which(patchxys$x%in%xset),which(patchxys$y%in%yset)) patchIDs<-unique(patchIDs) ctrl1$patches<-patchIDs # here ctrl2 is a one-off pesticide spray of patches that didn't have habitat removed #ctrl2<-control #ctrl2$stage<-ctrlStages[2] #ctrl2$type<-ctrlTypes[1] #ctrl2$intensity<-0.70 #pesticide has an average efficacy/percent kill of 70% of population #ctrl2$times[4]<-1 #ctrl2$patches<-1:numPatch #ctrl2$patches<-ctrl2$patches[-ctrl1$patches] for (repl in 1:numReps) { # set constant landscape where all cells are habitat K<-rep(meanKperCell,times=numPatch) # set landscape with random realisation of proportion in habitat and habitat fragmentation #K<-randomPatch(Lx=nx, Ly=ny, p = fragHab, A = propHab, dirs = 8, plt = FALSE) # allocate patches to control - only need to do here if random for each replicate # randomly sample same number of patches for removal as in swath simuations #nPatchCtrl<-round((1-(1/4))*numPatch) #ctrl1$patches<-sample(1:numPatch,nPatchCtrl,replace=FALSE) # This applies pesiticide to all patches that weren't randomly removed in ctrl1 #ctrl2$patches<-1:numPatch #ctrl2$patches<-ctrl2$patches[-ctrl1$patches] # initialise population (starting with egg masses) Egg<-rep(0,times=numPatch) # put initial occupied patch in middle of grid (lower left of central 4 squares if nx,ny are even integers) occpatch<-nx*(ny%/%2)+ceiling(nx/2) Egg[occpatch]<-rpois(initNumPatchOcc,numEggMasses*eggsPerFem) # write initial population to results array metaResArr[1,"sumFe",repl]<-numEggMasses metaResArr[1,"sumEggs",repl]<-sum(Egg) metaResArr[1,"numOccEggs",repl]<-length(which(Egg>0)) # time loop (note no zero-based arrays in R so have to bodge times) for (t in 2:(tMax+1)) { #ctrl type 1 kills all eggs and sets carrying capcity (K) to zero - non habitat if ((ctrl1$times[t-1]==1) & (ctrl1$stage=="E")) { Egg[ctrl1$patches]<-0 K[ctrl1$patches]<-0 } # eggs hatch into 1st instar larvae Templ<-numSurv(Egg,K,EtoLsurv) Larv<-rep(0,numPatch) # 1st instar larvae disperse through ballooning # loop through all patches and redistribute larvae according to probabilities of dispersal # w.r.t. distance between patches, actual number moved drawn from a multinomial distribution for (j in 1:numPatch) { # wrapped boundaries implemented by calculating distance between patches at grid edges as adjacent pdisp<-probDisp(Dists=distBWpatches[j,],dispScal=LdispDistScal,maxDist = maxLdispRad) Nemig<-rmultinom(1,Templ[j],pdisp) # if using absorbing boundaries add another bin for prob of dispersing off grid then draw from multinomial distribution and discard those off grid # pdisp2<-c(pdisp,1-sum(pdisp)) # Nemig<-rmultinom(1,Templ[j],pdisp2) Larv<-Larv+Nemig } # kill the larvae that land on non-habitat (K=0) patches Larv<-ifelse(K==0,0,Larv) # control type 2 applies pesticide to larvae in specific patches/grid cells # if ((ctrl2$times[t-1]==1) & (ctrl2$stage=="L")) { # Larv[ctrl2$patches]<-numSurv(Larv[ctrl2$patches],K[ctrl2$patches],(1-ctrl2$intensity)) # } # adults emerge and are allocated to male or female adults<-numSurv(Larv,K,LtoAsurv,dd=LtoAddcoef) Fe<-rbinom(length(adults),adults,sexRatio) M<-adults-Fe # females mate and lay eggs for (j in 1:numPatch) Fe[j]<-rbinom(1,Fe[j],probMated(M, Dists=distBWpatches[j,],distScal=MattrDistcoef,maxProb=MattrMult)) Egg<-numProd(Fe,K,eggsPerFem,dd=0) #write to results array metaResArr[t,"sumFe",repl]<-sum(Fe) metaResArr[t,"sumEggs",repl]<-sum(Egg) metaResArr[t,"numOccEggs",repl]<-length(which(Egg>0)) } # end time loop } #end replicate loop # write results array to file fname<-paste(numEggMasses,"em",swathwidth,"swRemHab",".csv",sep="") write.csv(as.data.frame.table(metaResArr, responseName = "value"), file = fname,row.names = FALSE) } # end sw loop #------------------------------------------------------------------------------