# Calculates by simulation the probability that a population # undergoing density-dependent, theta-logistic growth will fall # below the quasi-extinction threshold Nx at or before time tmax # Assumes that a population starts at its carrying capacity (K) ############# Simulation input parameters ########################## r <- 0.5967 # intrinsic rate of increase theta <- 1 # nonlinearity in density dependence sigma2 <- 0.0745 # environmental variance of r maxN <- 30000 # maximum initial population size (Nc) to simulate interval <- 100 # sample possible Nc values how often probcast <- 0.2 # per year probability of catastrophe alpha <- 0.5 # the shape1 (alpha) parameter from a Beta distribution characterizing the distribution of fish-kill severities beta <- 2.8 # the shape2 (beta) parameter from a Beta distribution characterizing the distribution of fish-kill severities Nx <- 50 # quasi-extinction threshold tmax <- 100 # Number of years being simulated NumReps <- 10000 # number of replicate population trajectories per Nc value ############# Population Growth Simulation ######################### sigma <- sqrt(sigma2) inter <- interval*(c(1:(maxN/interval))) M <- matrix(NA,nrow=length(inter),ncol=3) colnames(M) <- c("K","Prob.persist","Prob.extinct") for (j in 1:length(inter)) { N <- matrix(NA,nrow=tmax,ncol=NumReps) K <- j*interval Nc <- K N[1,] <- Nc for (i in 1:NumReps) { for (t in 2:tmax) { catastrophe <- 1 - (rbinom(1*1, size=1, prob=probcast)*rbeta(1*1, shape1=alpha, shape2=beta)) N[t,i] <- N[(t-1),i]*exp(r*(1-((N[(t-1),i]/K)^theta))+(sigma*rnorm(1,mean=0,sd=1)))*catastrophe } } M[j,1] <- K M[j,2] <- 1-((NumReps-length(which(N[tmax,]>0)))/NumReps) M[j,3] <- 1-M[j,2] } ############# Estimate of MVP ###################################### MVP <- interval*min(which(M[,2]>0.94999)) MVP write.table(M,"MVP.txt") ############# Plot of extinction probability ####################### plot(M[,1],M[,3], 'l', col = 'blue', ylim=c(0,1), xlab='Initial population size', ylab='Probability of extinction')