Exploring Spatial Patterns and Coexistance

Implementation of the Coupled Map Lattice model presented in “Self-Organization in Complex Ecosystems” SolĂ© & Bascompte (2006) Chp 3, page 87.

#load the libraries
library(plyr)
library(reshape2)
library(ggplot2)
library(viridis)

#Lotka-Voltera discrete time competition function
#@Ns are the population abundance
#@rs are the growth rate
#@as is the competition matrix
discrete<-function(Ns,rs,as){
  #return population abundance at time t+1
  return(as.numeric(Ns*exp(rs*(1-as%*%Ns))))
}

#an helper function returning line numbers of all neighbouring cells
find_neighbour<-function(vec,dat=dat){
  lnb<-which(dat[,"x"]%in%(vec["x"]+c(-1,-1,-1,0,1,1,1,0)) & dat[,"y"]%in%(vec["y"]+c(-1,0,1,1,1,0,-1,-1)))
  #remove own line
  lnb<-lnb[which(lnb!=vec["lnb"])]
  return(lnb)
}

Look at the dynamic for two species with equal growth rate, equal dispersion rate and equal interspecific competition.

set.seed(20160611)
generation<-1:60
dat<-expand.grid(x=1:50,y=1:50)
dat$N1<-0.5+rnorm(2500,0,0.1)
dat$N2<-0.5+rnorm(2500,0,0.1)
dat$lnb<-1:2500
list_dat<-list(dat)
#competition parameters
as<-matrix(c(1,1.2,1.2,1),ncol=2,byrow=TRUE)
#growth rate
rs<-c(1.5,1.5)
#dispersal rate
ds<-c(0.05,0.05)
#infos on neighbouring cells
neigh<-apply(dat,1,function(x) find_neighbour(x,dat))
#model the dynamic assuming absorbing boundary (ie ind that go out of the grid are lost to the system)
for(i in generation){  
  list_dat[[(i+1)]]<-rbind.fill(apply(list_dat[[i]],1,function(x){
    #population dynamics within one grid cell
    intern<-(1-ds)*discrete(Ns=x[c("N1","N2")],rs=rs,as=as)
    #grab the neighbouring cells
    neigh_cell<-list_dat[[i]][neigh[[x["lnb"]]],]
    #the number of immigrant coming into the focal grid cell
    imm<-(ds/8)*rowSums(apply(neigh_cell,1,function(y) discrete(Ns=y[c("N1","N2")],rs=c(1.5,1.5),as=as)))
    out<-data.frame(x=x["x"],y=x["y"],N1=intern[1]+imm[1],N2=intern[2]+imm[2],lnb=x["lnb"])
    return(out)
  }))
  #print(i)
}

#look at coexistence between the two species within one cell
cell525<-ldply(list_dat,function(x) c(x[525,"N1"],x[525,"N2"]))
cell525$Gen<-1:61
cell525<-melt(cell525,id.vars="Gen",measure.vars=1:2,variable.name="Species",value.name="Pop")
ggplot(cell525,aes(x=Gen,y=Pop,color=Species))+geom_path()

plot of chunk unnamed-chunk-2

#look at coexistence across the system
all<-ldply(list_dat,function(x) c(sum(x[,"N1"]),sum(x[,"N2"])))
all$Gen<-1:61
all<-melt(all,id.vars="Gen",measure.vars=1:2,variable.name="Species",value.name="Pop")
ggplot(all,aes(x=Gen,y=Pop,color=Species))+geom_path()

plot of chunk unnamed-chunk-2

#compare species-species abundance correlation at the beginning and at the end of the simulation
exmpl<-rbind(list_dat[[1]],list_dat[[60]])
exmpl$Gen<-rep(c("Gen : 1","Gen : 60"),each=2500)
ggplot(exmpl,aes(x=N1,y=N2))+geom_point()+facet_grid(.~Gen)

plot of chunk unnamed-chunk-2

#look at variation in spatial patterns at the beginning and the end of the simulation
exmplm<-melt(exmpl,id.vars=c("x","y","Gen"),measure.vars=c("N1","N2"),variable.name="Species",value.name="Pop")
ggplot(exmplm,aes(x,y,fill=Pop))+geom_raster()+scale_fill_viridis(option="plasma",direction=-1)+facet_grid(Gen~Species)

plot of chunk unnamed-chunk-2