Metapopulation model in 2 or 3D, using R / igraph

3d-metapop.r

# Metapop model in R

library(igraph)
library(plyr)
library(reshape2)
library(RColorBrewer)

colorize <- function(x, pal='Oranges')
{
   x <- x - min(x)
   x <- x / max(x)
   x <- round(x*100,0)
   cols <- colorRampPalette(brewer.pal(9, pal))(100)[x]
   return(cols)
}

createGeometricGraph <- function(n_sites, m_distance, dimensions)
{
   xy <- replicate(dimensions, runif(n_sites, 0, 1))
   A <- matrix(0, ncol=n_sites, nrow=n_sites)
   distance_matrix <- as.matrix(dist(xy, upper=T, diag=T))
   A[distance_matrix <= m_distance]  <- 1
   diag(A) <- 0
   return(list(xy = xy, graph = graph.adjacency(A)))
}

simulateMetaPop <- function(landscape, initial_occupancy, steps, migration_rate, extinction_rate, update_fraction)
{
   with(as.list(landscape),
        {
           ## Receivers
           ## Initiate population
           initial_state <- rbinom(length(V(graph)), 1, initial_occupancy)
           graph <- set.vertex.attribute(graph, 'state', V(graph), initial_state)
           graphs = list(graph)
           ## Start the loop
           for(now in c(1:steps))
           {
              ## Nodes to update
              n_updates <- round(update_fraction * length(V(graph)), 0)
              updated_nodes <- sample(V(graph), n_updates, replace=FALSE)
              for(patch in updated_nodes)
              {
                 neighbors <- graph[[patch]][[1]]
                 ## First we test for migration
                 if(runif(1) < migration_rate)
                 {
                    receiving_patch <- sample(neighbors, 1)
                    graph <- set.vertex.attribute(graph, 'state', receiving_patch, 1)
                 }
                 ## Then we test for extinction
                 if(runif(1) < extinction_rate) graph <- set.vertex.attribute(graph, 'state', patch, 0)
              }
              graphs[[now+1]] = graph
           }
           return(graphs)
        })
}



## 2D landscape
landscape = createGeometricGraph(80, 0.2, 2)

#plot(landscape$graph, layout=as.matrix(landscape$xy), vertex.label=NA, vertex.size=5, vertex.color='white', edge.arrow.size=0.2)
output <- simulateMetaPop(landscape, 0.1, 500, 0.2, 0.1, 0.3)
time_series_occupancy <- laply(output, function(x) get.vertex.attribute(x, 'state'))
average_occupancy <- aaply(time_series_occupancy, 2, mean)
plot(output[[20]], layout=as.matrix(landscape$xy), vertex.label=NA, vertex.size=5, vertex.color=colorize(sqrt(1-average_occupancy)), edge.arrow.size=0.2)


## 3D landscape
landscape = createGeometricGraph(80, 0.35, 3)
output <- simulateMetaPop(landscape, 0.1, 100, 0.2, 0.1, 0.3)
time_series_occupancy <- laply(output, function(x) get.vertex.attribute(x, 'state'))
average_occupancy <- aaply(time_series_occupancy, 2, mean)
rglplot(landscape$graph, layout=landscape$xy, vertex.color=colorize(sqrt(1-average_occupancy)), vertex.label=NA, vertex.size=7)