############################################################################### # Author: Che Elkin, 3 Nov. 2015 ############################################################################### # BIOL 410 # Metapopulation model # Incidence function model # Number of patches in system n.patches <- 27 # Patch X coordinate x.crd <- c(3.0527887, 6.3875999, 1.2157460, 5.8987072, 9.1664984, 4.2138941, 1.2712112, 1.7780712 ,8.1976360 ,1.2286834, 0.7647848, 4.2304007, 5.4847552, 7.2285807, 5.4960552, 0.2819228, 7.0801458, 1.0786554, 3.9252271, 9.7191894, 9.0732868, 7.2339198, 6.3941597, 4.9361327 ,1.6454546, 6.0316920 ,5.9592236) # Patch Y coordinate y.crd <- c(2.9720071, 10.7832916 , 3.4395872 , 4.5706549, 0.2968903 , 4.4312650, 6.0060003 , 9.3643911, 9.4566681 , 3.1582247, 5.6804270 , 9.2900524, 9.2236570 ,11.0599083 ,10.2841067, 5.9291035 , 8.5582610 ,0.9487339, 11.0514371, 7.1752323 , 2.7891322 , 7.5002027 , 7.9352540 , 4.6099435, 5.5533129, 11.8783566 , 5.9211253) # Patch area (m2) A <- c(7.06879832, 7.25512299, 1.08215340, 1.13369256, 2.17175274 ,1.12248860, 0.11057337, 0.32756816 ,0.46582475 ,0.08605073 ,0.37390697, 0.25872255, 1.83825572, 0.48647706, 5.74163538 ,0.53611531, 0.98189778 ,1.91684851, 2.12079187, 1.41999077, 0.31529620 ,2.93014128, 0.15878960, 2.91975299, 0.06806428, 0.44990127, 4.42327687) # Patch occupancy (2010) p <- c(1 ,1, 1, 0, 0, 1 ,1, 1 ,1 ,0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0,1, 0) # Model parameters alpha <- 1 # Exponential dispersal kernel parameter x <- 0.41 # Extinction parameter e <- 0.02 # Base extinction rate y <- 15.29 # Colonization parameter # Model functions d <- dist(cbind(x.crd,y.crd)) # Create distance matrix between patches edis <- as.matrix(exp(-alpha*d)) # Create realized connectivity based on disperal kernel diag(edis) <- 0 # Set the diagonal of the matrix to 0 E <- pmin(e/A^x,1) # Area based extinction rate # Plot dispersal kernel d.test <- c(1:12) plot(d.test,exp(-alpha*d.test), type="b",xlab="Distance (m)", ylab="Dispersal probability",ylim=c(0,0.5)) # Plot extinction probabilty as function of patch size A.test <- seq(0.001,8,0.1) plot(A.test,pmin(e/A.test^x,1), type="b",xlab="Patch size (m2)", ylab="Extinction probability",ylim=c(0,0.5)) # Plot the patches and plot(x.crd,y.crd,asp=1,xlab="",ylab="",cex=sqrt(A*2),pch=21,bg = 4*p,col="black") text(x.crd,y.crd, labels = c(1:length(A)), pos = 4) # Incidence function model metastep <- function(p,edis,E,y){ p <- p > 0 if(any(p)) { S <- rowSums(edis[,p,drop=FALSE]) C <- S^2/(S^2 + y) cond <- ifelse(p,(1-C)*E,C) p <- ifelse(runif(length(p)) < cond, !p, p) } } # Run the model for a set number of years tmax <- 200 # Number of years to simulate metapopulation dynamics occup <- matrix(0,nrow=length(p),ncol=tmax+1) occup[,1] <- p for(t in 1:tmax) { occup[, t +1 ] <- metastep(occup[,t], edis, E, y) } plot(colSums(occup), type = "l", col = "blue", lwd = 2, xlab = "Time", ylab = "Patches occupied") abline(h = mean(colSums(occup)), col = "red") # Plot the occupancy state of patches through time plot(x.crd,y.crd,asp=1,xlab="",ylab="",cex=sqrt(A*2),pch=21,bg = 4*occup[,1],col="black") for(t in 1:tmax) { points(x.crd,y.crd,asp=1,xlab="",ylab="",cex=sqrt(A*2),pch=21,bg = 4*occup[,t],col="black") Sys.sleep(0.01) }