#Models for Time to Extinction in Deteriorating Environments #Authors: Katherine Zarada and John M. Drake #Odum School of Ecology, University of Georgia #Parameters for the three models with two life histories and linear analysis #Model 1: Decreasing Birth rate (b1 varies) (b1 = 0.06for single trajectory plot in Fig 1) #Owl n=500: b0= 2, d0= 0.1, d1= .0038 #Owl n=2000: b0= 2, d0= 0.1, d1= 0.00095 #Rabbit n=500: b0= 12.5, d0= 0.5, d1= 0.024 #Rabbit n= 2000: b0= 12.5, d0= 0.5, d1= 0.006 #Linear decline n = 500: b0= 1, d0=0.8, d1= 0.00045 (b1 = 0.004 for the single trajectory plot in Fig 4) #Model 2: Increasing Death rate (d1 varies) (d1 = 0.6 for single trajectory plot in Fig 1) #Owl n=500: b0=2, d0= 0.1, b1= 0.0038 #Owl n=2000: b0= 2, d0=0.1, b1= 0.00095 #Rabbit n=500: b0= 12.5, d0= 0.5, b1= 0.024 #Rabbit n=2000: b0= 12.5, d0= 0.5, b1= 0.006 #Linear decline n = 500: b0= 1, d0=0.8, b1= 0.00045 (d1 = 0.004 for the single trajectory plot in Fig 4) #Model 3: Both rates are changing with time (b1=d1= varies) (b1 = 0.03 for single trajectory plot in Fig 1) #Owl n=500: b0=2, d0= 0.1, d2= 0.0038 #Owl n=2000: b0= 2, d0= 0.1, d2= 0.00095 #Rabbit n=500: b0=12.5, d0= 0.5, d2= 0.024 #Rabbit n= 2000: b0= 12.5, d0=0.5, d2 = 0.006 #Linear decline n = 500: b0= 1, d0=0.8, d2= 0.00045 (b1 = 0.002 for the single trajectory plot in Fig 4) #Birth rate model install.packages("adaptivetau") library(adaptivetau) parms <- list(b0=2, b1= 0.06 , d0= 0.1 , d1= 0.0095) #rabbit initial birth and death rate for b0 and d0, rate of deterioration for b1, and carrying capcity is determined through d1 transitions <- list(c(n=1), c(n=-1)) #n=1 is birth, n-1 is death f <- function(n, parms, t){ #birthrate model with birth rate decreasing exponentially with time birth <- n*parms$b0*exp(-1*parms$b1*t) death <- n*(parms$d0+parms$d1*n) return(c(birth, death)) #output is the birth and death rate, this changes as the population size (n) changes } set.seed(1) out <- ssa.adaptivetau(c(n=500), transitions, f, parms, tf=5000) sims <- 1000 #step is repeated 1000 times extinction.time <- c() #takes the formula above to get a probability for each transition at certain points, then gives an output of the time at which the population size is zero. for(i in 1:sims){ out <- ssa.adaptivetau(c(n=500), transitions, f, parms, tf=5000) extinction.time <- c(extinction.time, as.numeric(out[min(which(out[,2]==0)),1])) } hist(extinction.time, main='', xlab= 'Time to Extinction') #the output is a histogram of time to extinction birthrate<- summary(extinction.time) #gives mean for the distribution of time to extinctions sd(extinction.time) birthrate #Death rate model library(adaptivetau) parms <- list(b0=2, b1=.0038, d0= .1 , d1= 0.06) transitions <- list(c(n=1), c(n=-1)) f <- function(n, parms, t){ if(n<500) { birth <- n*(parms$b0-parms$b1*n) } else{ birth <- 0 } death <- n*(parms$d0*exp(parms$d1*t)) return(c(birth, death)) } set.seed(1) out <- ssa.adaptivetau(c(n=500), transitions, f, parms, tf=5000) sims <- 1000 extinction.time <- c() for(i in 1:sims){ out <- ssa.adaptivetau(c(n=2000), transitions, f, parms, tf=5000) extinction.time <- c(extinction.time, as.numeric(out[min(which(out[,2]==0)),1])) } hist(extinction.time, main='', xlab= 'Time to Extinction') deathrate<- summary(extinction.time) sd(extinction.time) deathrate #Both rates model library(adaptivetau) parms <- list(b0=2, b1=0.03, d0= .1, d2= 0.0038) transitions <- list(c(n=1), c(n=-1)) f <- function(n, parms, t){ birth <- n*(parms$b0*exp(-1*parms$b1*t)) death <- n*((parms$d0*exp(parms$b1*t))+ parms$d2*n) return(c(birth, death)) } set.seed(1) out <- ssa.adaptivetau(c(n=2000), transitions, f, parms, tf=5000) sims <- 1000 extinction.time <- c() for(i in 1:sims){ out <- ssa.adaptivetau(c(n=2000), transitions, f, parms, tf=5000) extinction.time <- c(extinction.time, as.numeric(out[min(which(out[,2]==0)),1])) } hist(extinction.time, main='', xlab= 'Time to Extinction') bothrate<- summary(extinction.time) sd(extinction.time) bothrate #Figures ############# #Fig 1 ############# #Figure 1 shows the deterministic carrying capacity decline for the three models # and a single trajectory of population decline for the three models. The code below # creates the deterministic carrying capacity by manipulating equations 2-4 in the paper # the code below those equations is the single trajectory code. The output is the population size at # various times until the time at which the population is extinct. The single trajectory code is similar to the model code above # but the output is the population size at time instead of the time at which the population is at zero. #Deterministic Carrying Capacity decline for Fig 1 t <- seq(from= 0, to= 60, by= 1) n <- ((2*exp(-0.06*t)) - .1) / 0.0038 # n = b0*exp(-b1*t)- d0/d1 n.2 <- (2 - (.1*exp(.06*t))) / 0.0038 #n = b0 - (d0*exp(d1*t))/b1 n.3 <- ((2*exp(-0.03*t)) - (.1*exp(0.03*t)))/0.0038 # n = (b0*exp(-b1*t) - d0*exp(b1*t))/d1 #single trajectories for Figure 1 #birth rate model, owl, 500, for 50 years to bifurcation library(adaptivetau) parms <- list(b0=2 , b1=.06 , d0= .1 , d1= 0.0038) #owl, 500 transitions <- list(c(n=1), c(n=-1)) f <- function(n, parms, t){ birth <- n*parms$b0*exp(-1*parms$b1*t) death <- n*(parms$d0+parms$d1*n) return(c(birth, death)) } set.seed(1) out.1 <- ssa.adaptivetau(c(n=500), transitions, f, parms, tf=5000) #Death rate model for same thing library(adaptivetau) parms <- list(b0=2, b1=.0038, d0= .1 , d1= 0.06) transitions <- list(c(n=1), c(n=-1)) f <- function(n, parms, t){ if(n<500) { birth <- n*(parms$b0-parms$b1*n) } else{ birth <- 0 } death <- n*(parms$d0*exp(parms$d1*t)) return(c(birth, death)) } set.seed(1) out.2 <- ssa.adaptivetau(c(n=500), transitions, f, parms, tf=5000) #both rates model library(adaptivetau) parms <- list(b0=2, b1= .03 , d0= .1 , d2= 0.0038) transitions <- list(c(n=1), c(n=-1)) f <- function(n, parms, t){ birth <- n*(parms$b0*exp(-1*parms$b1*t)) death <- n*((parms$d0*exp(parms$b1*t))+ parms$d2*n) return(c(birth, death)) } set.seed(1) out.3 <- ssa.adaptivetau(c(n=500), transitions, f, parms, tf=5000) #Figure 1 quartz() layout(matrix(c(1,1,2,2), 2,2, byrow=TRUE)) par(mar=c(4,4,1,1)) plot(out.1, type= 'l', xlab= "Time", ylab= 'Population Size', col= 'brown3', xlim= c(0,80)) abline(v=50, lty=2) lines(out.2, col= 'steelblue2') lines(out.3, col= 'black') mtext("a", side=3, line= -1, adj=.95, cex=0.9) plot(t, n, lwd= 3, type= 'l', xlim= c(0,80), xlab = 'Time', ylab = 'Carrying Capacity', col= 'brown3') lines(t, n.2, lwd=3, col= 'steelblue2') lines(t, n.3, lwd=3, col= 'black') abline(v=50, lty= 2) mtext("b", side=3, line= -1, adj=.95, cex=0.9) legend(60,400, c("Birth", "Death", "Both"), lty= c(1,1,1), lwd= c(2.5, 2.5,2.5), col= c("indianred3", "steelblue2", "black")) ############### # Fig 2 ############### #Figure 2 compares the delay for the different models at increasing levels of deterioration # the data below was obtained through running the models to get the mean and then subtracting the mean time # to extinction from the expected time to extinction #required package library(Hmisc) #data x <- c(0.001, 0.0015, 0.002, 0.003, 0.004, 0.005, 0.006, 0.011, 0.016, 0.021, 0.026, 0.031, 0.036, 0.041, 0.046, 0.051, 0.056, 0.061, 0.066, 0.071, 0.076, 0.081, 0.086, 0.091, 0.096, 0.101, 0.106, 0.111, 0.116, 0.121, 0.126, 0.131, 0.136, 0.141, 0.146, 0.151, 0.156, 0.161, 0.166, 0.171, 0.176, 0.181, 0.186, 0.191, 0.196, 0.201) W <- c(3.209711825, 2.463414297, 1.953971366, 1.506509077, 1.274397897, 1.113437964, 1.023945506, 0.730169912, 0.619173966, 0.557509551, 0.519245992, 0.493631543, 0.468965777, 0.457897805, 0.458530261, 0.446513606, 0.440505278, 0.406352679, 0.430069762, 0.432283356, 0.408882501, 0.435129406, 0.403190402, 0.405087768, 0.438924139, 0.414890829, 0.384849191, 0.407301363, 0.396233391, 0.408882501, 0.078740714, 0.378524636, 0.40761759, 0.409831185, 0.382319369, 0.411096096, 0.402241718, 0.391173747, 0.369986486, 0.392438658, 0.412993462, 0.413625918, 0.378208408, 0.394019796, 0.399395668, 0.386114102) SE.2 <- W*2 y.2 <- c(-1.732273554, 47.84515096, 66.13386322, 77.42257548, 80.66693161, 80.85354529, 80.21128774, 74.6607024, 70.7667329, 65.74608221, 63.37952794, 61.16347505, 59.48521462, 58.03335918, 57.97538536, 57.2601515, 56.00478083, 54.68963486, 55.11011707, 54.3365877, 53.60247008, 53.44565094, 52.9659038, 52.28986513, 52.80445548, 53.09928442, 51.19837478, 52.17142096, 51.27472178, 51.46188204, 51.45434704, 50.72181471, 51.04255681, 51.08367182, 50.3312858, 50.64071342, 51.11658799, 49.80296725, 50.03342004, 49.91109782, 51.2087939, 50.12899296, 49.93391251, 49.46553783, 50.44565167, 50.09585934) yplus.2 <-c(y.2+SE.2) yminus.2 <- c(y.2-SE.2) y.6 <- c(1.267726446, 11.84515096, 15.13386322, 17.42257548, 16.86693161, 16.95354529, 16.71128774, 14.7607024, 13.0667329, 11.94608221, 10.97952794, 10.36347505, 9.685214624, 9.243359182, 8.765385358, 8.430151499, 8.104780829, 7.85963486, 7.590117067, 7.346587696, 7.032470085, 6.885650944, 6.715903796, 6.509865126, 6.334455484, 6.20928442, 6.098374778, 5.941420959, 5.85472178, 5.721882037, 5.644347035, 5.501814706, 5.422556812, 5.373671819, 5.291285798, 5.23071342, 5.11658799, 5.102967245, 5.023420039, 4.951097815, 4.9187939, 4.848992964, 4.813912508, 4.765537835, 4.715651666, 4.705859336) SE.6 <- c(6.41942365, 0.987263086, 0.836106213, 0.631823077, 0.49963987, 0.428172395, 0.366191753, 0.239068191, 0.181514738, 0.151156872, 0.128388473, 0.116371818, 0.10751744, 0.094235874, 0.085381497, 0.079689397, 0.075262208, 0.07083502, 0.068937653, 0.064510464, 0.058818364, 0.057553453, 0.055656087, 0.051861354, 0.049963987, 0.048699076, 0.046801709, 0.045536798, 0.043639432, 0.041742065, 0.043006976, 0.039844699, 0.036682421, 0.037947332, 0.038579787, 0.038579787, 0.038579787, 0.039844699, 0.040477154, 0.04110961, 0.044904343, 0.046169254, 0.04806662, 0.050596443, 0.050596443, 0.056288542) yplus.6 <- c(y.6 + SE.6) yminus.6 <- c(y.6 - SE.6) y.10 <-c(30.13386322, 35.42257548, 35.76693161, 34.41128774, 33.33346581, 31.72677264, 30.45564387, 25.6303512, 22.38336645, 20.27304111, 18.65976397, 17.30173752, 16.25260731, 15.37667959, 14.59769268, 14.08007575, 13.31239041, 13.06481743, 12.52505853, 12.12329385, 11.81123504, 11.49782547, 11.1429519, 10.81993256, 10.52722774, 10.32964221, 10.11918739, 9.895710479, 9.66736089, 9.540941018, 9.332173518, 9.175907353, 8.936278406, 8.796835909, 8.650642899, 8.48035671, 8.358293995, 8.256483623, 8.096710019, 8.030548908, 7.88939695, 7.764496482, 7.686956254, 7.537768917, 7.457825833, 7.377929668) SE.10 <- c(1.706997481, 1.420495125, 1.066952483, 0.820294825, 0.652694109, 0.569842434, 0.505964426, 0.34532072, 0.258041857, 0.223256803, 0.18973666, 0.162541072, 0.153054239, 0.138507762, 0.128388473, 0.118269184, 0.118269184, 0.106884985, 0.101192885, 0.095500785, 0.096133241, 0.091073597, 0.08917623, 0.082219219, 0.076527119, 0.082219219, 0.079056942, 0.073997297, 0.068305197, 0.067672742, 0.069570109, 0.067040286, 0.061348187, 0.061348187, 0.058818364, 0.056920998, 0.057553453, 0.058818364, 0.056288542, 0.055656087, 0.05375872, 0.049963987, 0.049963987, 0.048699076, 0.049963987, 0.051228898) yplus.10 <- c(y.10 + SE.10) yminus.10 <- c(y.10 - SE.10) #Inset plot data i <- c(3750, 3529, 3333, 3158, 3000) #x values i.1 <- c(-34.67, -16.39, -9.59, -4.40, -1.73) #y values for birth i.2 <- c(-7.67, -5.39, -1.59, -0.40, 1.27)#death i.3 <- c(24.67, 25.8, 28.7, 28.3, 30.13)#both SE.1 <- c(16.41, 12.47, 15.46, 14.03, 6.42) SE.2 <- c(3.71, 3.12, 3.21, 2.98, 6.42) SE.3 <- c(4.97, 4.00, 3.91, 4.03, 1.71) i1.plus <- c(i.1+SE.1) i1.minus <- c(i.1-SE.1) i2.plus <- c(i.2 + SE.2) i2.minus <- c(i.2 - SE.2) i3.plus <- c(i.3 + SE.3) i3.minus <- c(i.3-SE.3) #plot code quartz() par(mar=c(5,4,1,2)+0.1) plot(x, y.2, ylim= c(0,90), ylab = "Delay", xlab= "Time to Transcritical Bifurcation", axes= F, col= 'red') errbar(x, y.2, yplus.2, yminus.2, col= 'brown3', add= TRUE) axis(side=1, at= c(0, 0.05, 0.1, 0.15, 0.2), label = c("3000", "60", "30", "20", "15")) axis(side=2, at= c(0, 20, 40, 60, 80) , label = c("0", "20", "40", "60", "80"), las=1) box() errbar(x, y.6, yplus.6, yminus.6, add= TRUE, col='steelblue2') errbar(x, y.10, yplus.10, yminus.10, add= TRUE) par(fig=c(0.4, 0.9, 0.6, 1), new=TRUE, mar=c(3.5,2,2,0.5)) #Code for inset plot(i, i.1, ylim= c(-30,10), xlim= c(3600, 3000), ylab= "", xlab= '') errbar(i, i.1, i1.plus, i1.minus, col= "brown3", add= TRUE) errbar(i, i.2, i2.plus, i2.minus, col= "steelblue2", add= TRUE) errbar(i, i.3, i3.plus, i3.minus, add= TRUE) abline(h=0, lty=2) ############ # Fig 3 ############ # Figure 3 compares the extinction delay for rabbit and owl birth rate models at a starting carrying capacity # of 500 or 2000 library(Hmisc) Q <- c(4.168830639, 3.17081581, 2.404595933, 1.890409585, 1.481527084, 1.415119253, 1.242458893, 0.889548706, 0.706136602, 0.640677454, 0.570158662, 0.546757807, 0.49963987, 0.484777165, 0.479085066, 0.458846488, 0.471811827, 0.460743855, 0.450940794, 0.437026773, 0.42184784, 0.396549619, 0.417420651, 0.420899157, 0.418053107, 0.407933818, 0.425958801, 0.38073823, 0.385797875, 0.414890829, 0.38073823, 0.405087768, 0.402874174, 0.384532963, 0.393387341, 0.384532963, 0.398763213, 0.364926842, 0.401293035, 0.382003141, 0.40761759, 0.376311042, 0.396549619, 0.367140436, 0.390225063, 0.404455313) SE <- Q*2 y <-c(-289.7322736, -146.154849, -87.86613678, -25.87742452, -1.233068388, 13.65354529, 24.11128774, 39.7607024, 40.8667329, 41.54608221, 42.47952794, 41.36347505, 40.68521462, 40.73335918, 40.17538536, 40.4801515, 39.91478083, 39.55963486, 39.49011707, 38.4565877, 38.63247008, 38.49565094, 38.3559038, 38.15986513, 38.13445548, 37.80928442, 38.49837478, 36.97142096, 37.14472178, 37.58188204, 37.14434704, 37.23181471, 37.53255681, 37.37367182, 37.0712858, 36.97071342, 36.59658799, 36.90296725, 36.63342004, 36.38109782, 36.8187939, 36.16899296, 36.47391251, 36.46553783, 36.26565167, 36.54585934) x <- c(0.001, 0.0015, 0.002, 0.003, 0.004, 0.005, 0.006, 0.011, 0.016, 0.021, 0.026, 0.031, 0.036, 0.041, 0.046, 0.051, 0.056, 0.061, 0.066, 0.071, 0.076, 0.081, 0.086, 0.091, 0.096, 0.101, 0.106, 0.111, 0.116, 0.121, 0.126, 0.131, 0.136, 0.141, 0.146, 0.151, 0.156, 0.161, 0.166, 0.171, 0.176, 0.181, 0.186, 0.191, 0.196, 0.201) yplus <- c(y+SE) yminus <- c(y-SE) #Delay for Owl n=2000 W <- c(3.209711825, 2.463414297, 1.953971366, 1.506509077, 1.274397897, 1.113437964, 1.023945506, 0.730169912, 0.619173966, 0.557509551, 0.519245992, 0.493631543, 0.468965777, 0.457897805, 0.458530261, 0.446513606, 0.440505278, 0.406352679, 0.430069762, 0.432283356, 0.408882501, 0.435129406, 0.403190402, 0.405087768, 0.438924139, 0.414890829, 0.384849191, 0.407301363, 0.396233391, 0.408882501, 0.078740714, 0.378524636, 0.40761759, 0.409831185, 0.382319369, 0.411096096, 0.402241718, 0.391173747, 0.369986486, 0.392438658, 0.412993462, 0.413625918, 0.378208408, 0.394019796, 0.399395668, 0.386114102) SE.2 <- 2*W y.2 <- c(-1.732273554, 47.84515096, 66.13386322, 77.42257548, 80.66693161, 80.85354529, 80.21128774, 74.6607024, 70.7667329, 65.74608221, 63.37952794, 61.16347505, 59.48521462, 58.03335918, 57.97538536, 57.2601515, 56.00478083, 54.68963486, 55.11011707, 54.3365877, 53.60247008, 53.44565094, 52.9659038, 52.28986513, 52.80445548, 53.09928442, 51.19837478, 52.17142096, 51.27472178, 51.46188204, 51.45434704, 50.72181471, 51.04255681, 51.08367182, 50.3312858, 50.64071342, 51.11658799, 49.80296725, 50.03342004, 49.91109782, 51.2087939, 50.12899296, 49.93391251, 49.46553783, 50.44565167, 50.09585934) yplus.2 <-c(y.2+SE.2) yminus.2 <- c(y.2-SE.2) #Delay for Rabbit n=500 E <- c(3.035470326, 2.137383471, 1.68454531, 1.273765442, 0.966392053, 0.85349874, 0.745348845, 0.475922788, 0.367772892, 0.301681289, 0.276699295, 0.251084846, 0.228316447, 0.20523182, 0.194480076, 0.181830965, 0.178668688, 0.170446766, 0.15431915, 0.153054239, 0.154002922, 0.147994594, 0.137875306, 0.132499434, 0.134080573, 0.126174879, 0.122063918, 0.122380145, 0.119850323, 0.115739362, 0.119217868, 0.108782352, 0.113525768, 0.11036349, 0.106884985, 0.108782352, 0.115423135, 0.111312174, 0.101192885, 0.109414807, 0.112260857, 0.098663063, 0.103406479, 0.098346835, 0.100876657, 0.098979291) SE.3 <- 2*E y.3 <- c(-537.8758249 ,-338.9172166, -302.4379124, -143.4586083, -97.81895622, -73.07516497, -55.07930414, -19.32507499, -7.879739054, -3.479801184, 0.597083659, 2.465295972, 3.996782643, 4.64083354, 4.96443859, 5.734787748, 6.240074556, 6.311543855, 6.379154169, 6.723720777, 6.876370725, 6.920792286, 7.111211339, 7.287738188, 7.190043491, 7.369942328, 7.123246935, 7.571118695, 7.501070475, 7.547720456, 7.493366469, 7.238428818, 7.561795405, 7.421093441, 7.402905309, 7.672941557, 7.78618061, 7.636982454, 7.409181778, 7.426164767, 7.630932813, 7.496155664, 7.574215995, 7.38724699, 7.427164159, 7.565692414) yplus.3 <- c(y.3+ SE.3) yminus.3 <- c(y.3 - SE.3) #Delay for Rabbit n=2000 P <- c(3.209711825, 1.625726945, 1.28388473, 0.942042515, 0.773176888, 0.670402864, 0.57616699, 0.376311042, 0.314330399, 0.248238796, 0.22926513, 0.207129187, 0.182147193, 0.176455093, 0.165070894, 0.155584061, 0.146729683, 0.152738011, 0.140405128, 0.132499434, 0.126491106, 0.119850323, 0.126174879, 0.125226195, 0.122063918, 0.117320501, 0.110995946, 0.11036349, 0.110047263, 0.107201213, 0.104671391, 0.104671391, 0.102141568, 0.102774024, 0.096765696, 0.101825341, 0.102457796, 0.097398152, 0.104671391, 0.096765696, 0.096765696, 0.096449469, 0.091389824, 0.09486833, 0.089808686, 0.096133241) SE.4 <- 2*P y.4 <- c(-224.8758249, -111.9172166, -70.43791243, -32.95860829, -16.81895622, -5.675164974, -2.179304145, 9.374925012, 12.92026095, 14.22019882, 14.29708366, 14.76529597, 14.18678264, 14.57083354, 14.21443859, 13.70478775, 13.71007456 ,13.84154385, 13.46915417, 13.26372078, 13.43637073, 12.84079229, 12.85121134, 13.05773819, 12.85004349, 12.68994233, 12.54324694, 12.48111869, 12.14107048, 12.17772046, 12.09336647, 12.13842882, 12.06179541, 12.00109344, 11.83290531, 11.66294156, 11.69618061, 11.56698245, 11.74918178, 11.47616477, 11.57093281, 11.36615566, 11.264216, 11.46724699, 11.29716416, 11.30569241) yplus.4 <- c(y.4+ SE.4) yminus.4 <- c(y.4 - SE.4) #delay combined for owl and rabbit birth rate decreasing for the n =2000 and n= 500 quartz() library(Hmisc) plot(x, y, type= 'l', ylim= c(0,85), ylab = "Delay", xlab= "Time to Transcritical Bifurcation", axes= F, col= "brown3") par(mar=c(4.5,4,1,1)) errbar(x, y, yplus, yminus, add= TRUE, col= "brown3") axis(side=1, at= c(0, 0.05, 0.1, 0.15, 0.2), label = c("3000", "60", "30", "20", "15")) axis(side=2, at= c(0, 20, 40, 60, 80) , label = c("0", "20", "40", "60", "80"), las=1) box() lines(x, y.2, col= 'steelblue2') errbar(x, y.2, yplus.2, yminus.2, add= TRUE, col='steelblue2') lines(x, y.3, pch= 17, bg= "brown3") errbar(x, y.3, yplus.3, yminus.3, add= TRUE, pch=17, col= 'brown3') lines(x, y.4, col= 'steelblue2') errbar(x, y.4, yplus.4, yminus.4, add= TRUE, col= 'steelblue2', pch= 17) legend(0.14, 80, legend= c("Owl", "Rabbit"), pch= c(16, 17), bty='n') ############ #Fig 4 ############ #Figure 4 is similar to Fig 1 and Fig 2 but is looking at the approximately linear analysis results #data for Fig 4a #birth rate rate = c(0.000233, 0.000234888, 0.000247937, 0.000262522, 0.000278929, 0.000297525, 0.000318777, 0.000343298, 0.000371906, 0.000405716, 0.000446287, 0.000495875, 0.000557859, 0.000637553, 0.000743812, 0.000892574, 0.000991749, 0.001115718, 0.001275106, 0.001487624, 0.001785148, 0.002231436, 0.002625218, 0.002975247, 0.003432978, 0.004057155, 0.004462871, 0.004958746, 0.00637553, 0.008925742, 0.014876237, 0.022314355) br = c(-114.9, -69.5, -58.7, -48, -43.3, -31.6, -23.6, -13.6, -7.8, 2.6, 7.2, 16.9, 20.6, 27, 34.6, 36, 38.7, 43.1, 41.9, 43.1, 43.5, 45.2, 42.1, 41.5, 41, 39.07, 37.94, 37.55, 35.35, 32.29, 27.61, 24.26) SE.b = c(2.598285439, 2.758454803, 2.61868213, 2.496301985, 2.264190805, 2.212013223, 2.067813362, 2.006781403, 1.875863108, 1.726287375, 1.732928158, 1.563430075, 1.424922314, 1.266808431, 1.142530919, 0.98663063, 0.93856401, 0.876583367, 0.790569415, 0.673881369, 0.610003361, 0.549603857, 0.481931115, 0.432599584, 0.389592608, 0.356388692, 0.332355382, 0.301048833, 0.255512035, 0.227367764, 0.161592388, 0.130602067) brplus = br+ SE.b brminus = br-SE.b #death rate dr = c(-105.3, -58.9, -52.2, -44.1, -35.2, -29.1, -21.9, -11.5, -5.7, -1.2, 6.7, 13.5, 18, 25, 29, 32.8, 34, 38.7, 36.6, 37.5, 37.7, 38, 35.8, 35.6, 34.32, 33.05, 32.25, 31.41, 29.12, 26.21, 21.58, 18.41) SE.d = c(2.386254722, 2.398587605, 2.29739472, 2.160468097, 2.064967312, 1.995080976, 1.880290297, 1.810087733, 1.671263743, 1.546353776, 1.465715695, 1.386974982, 1.270286936, 1.145693196, 1.000544652, 0.895557033, 0.804167209, 0.759579094, 0.669770408, 0.590397239, 0.522092042, 0.445248695, 0.39496848, 0.358602287, 0.320654955, 0.296937872, 0.266580007, 0.250136163, 0.208710326, 0.166652033, 0.115423135, 0.08633018) drplus = dr + SE.d drminus = dr-SE.d #both rates bor = c(-109.4, -68.5, -60.7, -47.2, -38.1, -32.8, -20.5, -16, -8, 0.7, 8.4, 12.7, 18.5, 25.9, 31.4, 35.1, 35.5, 39.9, 38.7, 39.6, 39.6, 40.9, 38.9, 38, 37.4, 36.18, 35.04, 34.29, 31.49, 29.1, 24.11, 20.84) SE.bo = c(2.512745829, 2.592435226, 2.463414297, 2.313838564, 2.180390447, 2.168057564, 2.000456848, 1.849616203, 1.743363674, 1.660511999, 1.546037548, 1.461288507, 1.362309216, 1.233920743, 1.082763871, 0.945837248, 0.880694328, 0.817448775, 0.7402892, 0.646053326, 0.556877096, 0.477820154, 0.413625918, 0.396233391, 0.364294386, 0.325714599, 0.302313744, 0.279229117, 0.230213814, 0.190052887, 0.138191534, 0.106568757) borplus = bor + SE.bo borminus = bor - SE.bo #single trajectory for Fig 4b t <- seq(from= 0, to= 60, by= 1) n.4 <- ((1*exp(-0.004*t)) - .8) / 0.00045 # n = b0*exp(-b1*t)- d0/d1 n.5 <- (1 - (0.8*exp(.004*t))) / 0.00045 #n = b0 - (d0*exp(d1*t))/b1 n.6 <- ((1*exp(-0.002*t)) - (.8*exp(0.002*t)))/0.00045 # n = (b0*exp(-b1*t) - d0*exp(b1*t))/d1 quartz() layout(matrix(c(1,1,2,2), 2,2, byrow=TRUE)) par(mar=c(4.5,4.5,1,1)) plot(rate, br,type ="l", lwd= 2, ylim= c(0,50), ylab = "Delay", xlab= "Time to Transcritical Bifurcation", axes= F, col= 'red') errbar(rate, br, brplus, brminus, col= 'brown3', add= TRUE) axis(side=1, at= c(0.000233, 0.00637553,0.014876237, 0.022314355, 0.04462871), label = c("1000","35","15", "10", "5")) axis(side=2, at= c(0, 20, 40, 60, 80) , label = c("0", "20", "40", "60", "80"), las=1) box() errbar(rate, dr, drplus, drminus, add= TRUE, col='steelblue2') lines(rate, dr, lwd= 2, col= 'steelblue2') errbar(rate, bor, borplus, borminus, add= TRUE) lines(rate, bor, lwd=2) mtext("a", side=3, line= -1, adj=.95, cex=0.9) plot(t, n.4, lwd= 3, type= 'l', xlim= c(0,55), ylim = c(-10,500), xlab = 'Time', ylab = 'Carrying Capacity', col= 'brown3') lines(t, n.5, lwd=3, col= 'steelblue2') lines(t, n.6, lwd=3, col= 'black') abline(v=50, lty= 2) mtext("b", side=3, line= -1, adj=.95, cex=0.9) legend(45,440, c("Birth", "Death", "Both"), lty= c(1,1,1), lwd= c(2.5, 2.5,2.5), col= c("indianred3", "steelblue2", "black")) ############ #Fig 5 ############ # Figure 5 is examining the extinction debt. The basic model has been changed to output the # number of individuals remaining at the deterministic time of extinction rates.2 = c(0.001, 0.0015, 0.002, 0.003, 0.004, 0.005, 0.006, 0.011, 0.016, 0.021, 0.026, 0.031, 0.036, 0.041, 0.046, 0.051, 0.056, 0.061, 0.066, 0.071, 0.076, 0.081, 0.086, 0.091, 0.096, 0.101, 0.106, 0.111, 0.116, 0.121, 0.126, 0.131, 0.136, 0.141, 0.146, 0.151, 0.156, 0.161, 0.166, 0.171, 0.176, 0.181, 0.186, 0.191, 0.196, 0.201) ext.times.1 = c(2995.732274, 1997.154849, 1497.866137, 998.5774245, 748.9330684, 599.1464547, 499.2887123, 272.3392976, 187.2332671, 142.6539178, 115.2204721, 96.63652495, 83.21478538, 73.06664082, 65.12461464, 58.7398485, 53.49521917, 49.11036514, 45.38988293, 42.1934123, 39.41752992, 36.98434906, 34.8340962, 32.92013487, 31.20554452, 29.66071558, 28.26162522, 26.98857904, 25.82527822, 24.75811796, 23.77565296, 22.86818529, 22.02744319, 21.24632818, 20.5187142, 19.83928658, 19.20341201, 18.60703275, 18.04657996, 17.51890218, 17.0212061, 16.55100704, 16.10608749, 15.68446217, 15.28434833, 14.90414066) n.vec.3 = rep(NA, length(rates.2)) for(i in 1:length(rates.2)){ parms <- list(b0=2 , b1=rates.2[i] , d0= .1 , d1= 0.00095) #owl, 500 transitions <- list(c(n=1), c(n=-1)) f <- function(n, parms, t){ birth <- n*parms$b0*exp(-1*parms$b1*t) death <- n*(parms$d0+parms$d1*n) return(c(birth, death)) } set.seed(1) out.1 <- ssa.adaptivetau(c(n=500), transitions, f, parms, tf=5000) time = out.1[,1] x= findInterval(ext.times.1[i], time) n= as.numeric(out.1[x[1],2]) n.vec.3[i] = n } n.vec.4 = rep(NA, length(rates.2)) for(i in 1:length(rates.2)){ parms <- list(b0=2, b1=0.00095, d0= .1, d1= rates.2[i]) transitions <- list(c(n=1), c(n=-1)) f <- function(n, parms, t){ if(n<500) { #need to change this to reflect the initial pop size birth <- n*(parms$b0-parms$b1*n) } else{ birth <- 0 } death <- n*(parms$d0*exp(parms$d1*t)) return(c(birth, death)) } set.seed(1) out.2 <- ssa.adaptivetau(c(n=500), transitions, f, parms, tf=5000) #also needs to be the initial pop size time = out.2[,1] x= findInterval(ext.times.1[i], time, all.inside = TRUE) n= as.numeric(out.2[x[1],2]) n.vec.4[i] = n } rates.3 = c(0.0005, 0.00075, 0.001 ,0.0015 ,0.002, 0.0025, 0.003, 0.0055, 0.008, 0.0105, 0.013, 0.0155 ,0.018, 0.0205, 0.023, 0.0255, 0.028, 0.0305, 0.033, 0.0355, 0.038, 0.0405, 0.043, 0.0455, 0.048, 0.0505, 0.053, 0.0555 ,0.058, 0.0605, 0.063, 0.0655, 0.068, 0.0705,0.073, 0.0755, 0.078, 0.0805, 0.083, 0.0855, 0.088, 0.0905, 0.093, 0.0955, 0.098, 0.1005) n.vec.5 = rep(NA, length(rates.3)) for(i in 1:length(rates.3)){ parms <- list(b0=2, b1=rates.3[i], d0= .1, d2= 0.00095) transitions <- list(c(n=1), c(n=-1)) f <- function(n, parms, t){ birth <- n*(parms$b0*exp(-1*parms$b1*t)) death <- n*((parms$d0*exp(parms$b1*t))+ parms$d2*n) return(c(birth, death)) } set.seed(1) out.3 <- ssa.adaptivetau(c(n=500), transitions, f, parms, tf=5000) time = out.3[,1] x= findInterval(ext.times.1[i], time, all.inside = TRUE) n= as.numeric(out.3[x[1],2]) n.vec.5[i] = n } #Plot for Fig 5 quartz() Birth = n.vec.3[1:7] Death = n.vec.4[1:7] Both = n.vec.5[1:7] debt = as.matrix(rbind(Birth, Both, Death)) barplot(debt, xlab = "Time to Transcritical Bifurcation", ylab = "Extinction Debt", col = c("brown3","black", "steelblue2"), names.arg = c("3000", "2000", "1500", "1000", "750", "600", "500"), axes = TRUE, beside= TRUE, axisnames = TRUE)