# Individual-based northern quoll population model reported in (Kelly & Phillips 2016) # load parameter files load("Kept_sims_poponlyhalfpcED.RData") pars<-apply(gold[,1:6], 2, mean) # parameters from ABC computation ### Functions # Initialises a matrix with n individuals init.inds<-function(K, prob.attack){ n<- K*2 S<-rbinom(n, 1, 0.5) #sex A<-rep(1, n) #age P<-rbinom(n, 1, prob.attack) # toad attack probability cbind(S, A, P) } # Density Dependance function ddep.fun<- function(x, K){ alpha<- exp(x*(log(0.3)/K)) # calculates shape of the curve alpha } # Offspring number function # input data maxSize<-8 # max litter size trueProb<-0.2 X<-c(8,5,7,7,8,7,2,5,2,8,6,8,8,7) # raw litter size data # fit binomial function fitBinomial <- function(X, maxSize) { binLL<-function(estP){-sum(dbinom(x=X, size=maxSize, prob = estP, log = TRUE))} optimize(f = binLL, lower=0, upper=1) } fitBinomial(X, maxSize) lit.prob<-fitBinomial(X, maxSize) lit.prob<-lit.prob$minimum # Mate Selection function mating<- function(female, male, maxmates){ if (nrow(female)==0) return(NULL) mates<-nrow(male)*maxmates if (mates>nrow(female)) {mates<-nrow(female)} mated<- female[sample(nrow(female),mates),] mated } # Function that reproduces individuals # density dependance acting on % FEMALE BREEDING according to K repro<-function(popmat, K, prob.attack){ if (length(popmat[,1])==0) return(popmat) # collect breeding females and find their mates male<-subset(popmat, popmat[,"S"]==1 & popmat[,"A"]>0) if (length(male)==0) return(popmat) female<-subset(popmat, popmat[,"S"]==0 & popmat[,"A"]>0) #matrix of all adult females female<-mating(female, male, 10) if (is.null(nrow(female))) return(popmat) # breeding females based on density dependance: ddep<- ddep.fun(x=nrow(female), K) fdens<- as.integer(nrow(female)*ddep) if (fdens==1) fdens<-2 female<-female[sample(nrow(female), fdens), ] off.no<- rbinom(nrow(female),8, lit.prob) # 8 max no. of offspring # get info for next generation and combine offspring<-matrix(female[rep(1:nrow(female), off.no),], ncol=ncol(female), dimnames=list(NULL, colnames(female))) if (is.null(nrow(offspring)) | is.na(nrow(offspring))) return(popmat) offspring[,"S"]<-rbinom(nrow(offspring), 1, 0.5) # random sex allocation offspring[,"A"]<-rep(0, nrow(offspring)) # age=0 offspring[,"P"]<-rbinom(nrow(offspring), 1, prob.attack) # offspring toad attack probability popmat<-rbind(popmat, offspring) popmat } # Male die off mdieoff<-function(popmatrix, msurv){ psurv<-rep(1, length(popmatrix[,1])) psurv[which(popmatrix[,"S"]==1 & popmatrix[,"A"]==1)]<-msurv #one year old males psurv[which(popmatrix[,"S"]==1 & popmatrix[,"A"]>1)]<-0 #two year old males popmatrix<-subset(popmatrix, rbinom(length(popmatrix[,1]), 1, psurv)==1) # probabilistic survival popmatrix } # Survival (based on age and selection) age<-function(popmatrix, selection=F, fsurv1, fsurv2){ if (length(popmatrix[,1])==0) return(popmatrix) # Survival of adults Ad<-subset(popmatrix, popmatrix[,"A"]>0) #Grabs all adults psurv<-rep(0, length(Ad[,1])) psurv[which(Ad[,"S"]==0 & Ad[,"A"]==1)]<-fsurv1 #one year old females psurv[which(Ad[,"S"]==0 & Ad[,"A"]==2)]<-fsurv2 #two year old females psurv[which(Ad[,"S"]==1 & Ad[,"A"]==1)]<-1 #one year old males (already gone through die off) Ad<-subset(Ad, rbinom(length(Ad[,1]), 1, psurv)==1) # probabilistic survival Juv<-subset(popmatrix, popmatrix[,"A"]==0) popmatrix<-rbind(Ad, Juv) # Survival if toads are present if (selection==T) popmatrix<-subset(popmatrix, popmatrix[,"P"]==0) # 0==does not attack = surviving quolls # Gather survivors, age them and return popmatrix[,"A"]<-popmatrix[,"A"]+1 popmatrix } ### Mother function # dem.pars is a vector with alpha, fs1, fs2, msurv, beta, prob.d mother<-function(dem.pars=pars, K=1000, prob.attack, gens, sel.time){ fsurv1<-dem.pars[2] fsurv2<-dem.pars[3] msurv<-dem.pars[4] pop<-init.inds(K, prob.attack=prob.attack) popsize<-nrow(pop) sel<-FALSE for (g in 2:gens){ pop<-repro(popmat=pop, K=K, prob.attack=prob.attack) # females reproduce (density dep) pop<-mdieoff(pop, msurv) # males die off if (g>sel.time) sel<-TRUE # have toads arrived? pop<-age(pop, selection=sel, fsurv1=fsurv1, fsurv2=fsurv2) # everyone gets a little older if (length(pop[,1])==0) { # did we go extinct? pe<-TRUE break } else {pe<-FALSE } popsize<-c(popsize, length(pop[,1])) # new population size appended to popsize vector } return(pe) } ### Run model probs<-seq(from = 0, to = 1, by = 0.01) pe<-vector() n<-100 for (j in probs){ for (i in 1:n){ prob.ex<-mother(K=1000, prob.attack=j, gens=100, sel.time=30) pe[i]<-prob.ex } print(length(pe[pe==TRUE])/length(pe)) }