**R – biologyforfun**, and kindly contributed to R-bloggers)

Today is a rainy day and I had to drop my plans for going out hiking, instead I continued reading “Self-Organization in Complex Ecosystems” from Richard Solé and Jordi Bascompte.

As I will be busy in the coming weeks with spatial models at the iDiv summer school I was closely reading chapter three on spatial self-organization and decided to try and implement in R one of the Coupled Map Lattice Models that is described in the book.

The model is based on a discrete time version of the Lotka-Volterra competition model, for example for two competing species we model the population normalized abundances (the population abundance divided by the carrying capacity) N1 and N2 using the following equations:

Where r is the growth rate and the are the interspecific competition coefficients. Note that since we use normalized abundance the intraspecific competition coefficient is equal to one.

Space can be introduced in this model by considering dispersion between neighboring cells. Say we have a 50 x 50 lattice map and that individuals disperse between neighboring cells with a probability D, as a result the normalized population abundances in cell k at time t+1 will be:

Where j is an index of all 8 neighboring cells around cell k. We have a similar equation for the competing species.

Now let’s implement this in R (you may a nicer version of the code here).

I start by defining two functions: one that will compute the population dynamics based on the first two equations, and the other one that will return for each cell in lattice the row numbers of the neighboring cells:

#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) }

Then I create a lattice, fill it with the two species and some random noise.

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)

The next step is to initialize all the parameters, in the example given in the book the two species have equal growth rate, equal dispersal and equal interspecific competition.

#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))

Now we are ready to start the model (the code is rather ugly, sorry for that, saturday coding makes me rather lazy …)

generation<-1:60 #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) }

First let’s reproduce figure 3.8 from the book:

#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() #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()

And now figure 3.9 plus an extra figure looking at the correlation in population abundance between the two species at the beginning and the end of the simulation:

#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) #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)

Neat stuff, I am looking forward to implement other models in this book!

Filed under: Biological Stuff, R and Stat Tagged: ecology, R

**leave a comment**for the author, please follow the link and comment on their blog:

**R – biologyforfun**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...