<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="http://www.r-bloggers.com/exploring-spatial-patterns-and-coexistance/"></div></div><div><a data-url="http://www.r-bloggers.com/exploring-spatial-patterns-and-coexistance/" data-counturl="http://www.r-bloggers.com/exploring-spatial-patterns-and-coexistance/" data-count="horizontal" data-text="Exploring Spatial Patterns and Coexistance" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://biologyforfun.wordpress.com/2016/06/11/exploring-spatial-patterns-and-coexistance/"> R – biologyforfun</a></strong>, and kindly contributed to <a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Today is a rainy day and I had to drop my plans for going out hiking, instead I continued reading <a target="_blank" href="http://press.princeton.edu/titles/8224.html">“Self-Organization in Complex Ecosystems”</a> from Richard Solé and Jordi Bascompte.</p>
<p>As I will be busy in the coming weeks with spatial models at the <a target="_blank" href="https://www.idiv.de/ydiv/summer_school_2016.html">iDiv summer school</a> 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.</p>
<p>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:</p>
<p><img title="N_{1,t+1} = \phi_{N_{1}}(N_{1,t},N_{2,t}) = N_{1,t} * exp[r_{N_{1}} * (1 - N_{1} - \alpha_{12} * N_{2})] " alt="N_{1,t+1} = \phi_{N_{1}}(N_{1,t},N_{2,t}) = N_{1,t} * exp[r_{N_{1}} * (1 - N_{1} - \alpha_{12} * N_{2})] " class="latex" src="https://s0.wp.com/latex.php?latex=N_%7B1%2Ct%2B1%7D+%3D+%5Cphi_%7BN_%7B1%7D%7D%28N_%7B1%2Ct%7D%2CN_%7B2%2Ct%7D%29+%3D+N_%7B1%2Ct%7D+%2A+exp%5Br_%7BN_%7B1%7D%7D+%2A+%281+-+N_%7B1%7D+-+%5Calpha_%7B12%7D+%2A+N_%7B2%7D%29%5D+&bg=f9f7f5&%23038;fg=444444&%23038;s=0"></p>
<p><img title="N_{2,t+1} = \phi_{N_{2}}(N_{1,t},N_{2,t}) = N_{2,t} * exp[r_{N_{2}} * (1 - N_{2} - \alpha_{21} * N_{1})] " alt="N_{2,t+1} = \phi_{N_{2}}(N_{1,t},N_{2,t}) = N_{2,t} * exp[r_{N_{2}} * (1 - N_{2} - \alpha_{21} * N_{1})] " class="latex" src="https://s0.wp.com/latex.php?latex=N_%7B2%2Ct%2B1%7D+%3D+%5Cphi_%7BN_%7B2%7D%7D%28N_%7B1%2Ct%7D%2CN_%7B2%2Ct%7D%29+%3D+N_%7B2%2Ct%7D+%2A+exp%5Br_%7BN_%7B2%7D%7D+%2A+%281+-+N_%7B2%7D+-+%5Calpha_%7B21%7D+%2A+N_%7B1%7D%29%5D+&bg=f9f7f5&%23038;fg=444444&%23038;s=0"></p>
<p>Where r is the growth rate and the <img title="\alpha" alt="\alpha" class="latex" src="https://s0.wp.com/latex.php?latex=%5Calpha&bg=f9f7f5&%23038;fg=444444&%23038;s=0"> are the interspecific competition coefficients. Note that since we use normalized abundance the intraspecific competition coefficient is equal to one.</p>
<p>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:</p>
<p><img title="N_{1,t+1,k} = (1 - D) * \phi_{N_{1}}(N_{1,t,k},N_{2,t,k}) + \frac{D}{8} * \sum_{j=1}^{8} \phi_{N_{1}}(N_{1,t,j},N_{2,t,j}) " alt="N_{1,t+1,k} = (1 - D) * \phi_{N_{1}}(N_{1,t,k},N_{2,t,k}) + \frac{D}{8} * \sum_{j=1}^{8} \phi_{N_{1}}(N_{1,t,j},N_{2,t,j}) " class="latex" src="https://s0.wp.com/latex.php?latex=N_%7B1%2Ct%2B1%2Ck%7D+%3D+%281+-+D%29+%2A+%5Cphi_%7BN_%7B1%7D%7D%28N_%7B1%2Ct%2Ck%7D%2CN_%7B2%2Ct%2Ck%7D%29+%2B+%5Cfrac%7BD%7D%7B8%7D+%2A+%5Csum_%7Bj%3D1%7D%5E%7B8%7D+%5Cphi_%7BN_%7B1%7D%7D%28N_%7B1%2Ct%2Cj%7D%2CN_%7B2%2Ct%2Cj%7D%29+&bg=f9f7f5&%23038;fg=444444&%23038;s=0"></p>
<p>Where j is an index of all 8 neighboring cells around cell k. We have a similar equation for the competing species.</p>
<p>Now let’s implement this in R (you may a nicer version of the code <a target="_blank" href="http://rpubs.com/hughes/188614">here</a>).</p>
<p>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:</p>
<pre>#Lotka-Voltera discrete time competition function
#@Ns are the population abundance
#@rs are the growth rate
#@as is the competition matrix
discrete&lt;-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&lt;-function(vec,dat=dat){
  lnb&lt;-which(dat[,&quot;x&quot;]%in%(vec[&quot;x&quot;]+c(-1,-1,-1,0,1,1,1,0)) &amp;
+ dat[,&quot;y&quot;]%in%(vec[&quot;y&quot;]+c(-1,0,1,1,1,0,-1,-1)))
  #remove own line
  lnb&lt;-lnb[which(lnb!=vec[&quot;lnb&quot;])]
  return(lnb)
}
</pre>
<p>Then I create a lattice, fill it with the two species and some random noise.</p>
<pre>dat&lt;-expand.grid(x=1:50,y=1:50)
dat$N1&lt;-0.5+rnorm(2500,0,0.1)
dat$N2&lt;-0.5+rnorm(2500,0,0.1)
dat$lnb&lt;-1:2500
list_dat&lt;-list(dat)</pre>
<p>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.</p>
<pre>#competition parameters
as&lt;-matrix(c(1,1.2,1.2,1),ncol=2,byrow=TRUE)
#growth rate
rs&lt;-c(1.5,1.5)
#dispersal rate
ds&lt;-c(0.05,0.05)
#infos on neighbouring cells
neigh&lt;-apply(dat,1,function(x) find_neighbour(x,dat))</pre>
<p>Now we are ready to start the model (the code is rather ugly, sorry for that, saturday coding makes me rather lazy …)</p>
<pre>generation&lt;-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)]]&lt;-rbind.fill(apply(list_dat[[i]],1,function(x){
    #population dynamics within one grid cell
    intern&lt;-(1-ds)*discrete(Ns=x,rs=rs,as=as)
    #grab the neighbouring cells
    neigh_cell&lt;-list_dat[[i]][neigh[[x[&quot;lnb&quot;]]],]
    #the number of immigrant coming into the focal grid cell
    imm&lt;-(ds/8)*rowSums(apply(neigh_cell,1,function(y) {
    discrete(Ns=y,rs=c(1.5,1.5),as=as)}))
    out&lt;-data.frame(x=x[&quot;x&quot;],y=x[&quot;y&quot;],N1=intern[1]+imm[1],N2=intern[2]
    + +imm[2],lnb=x[&quot;lnb&quot;])
    return(out)
  }))
  #print(i)
}</pre>
<p>First let’s reproduce figure 3.8 from the book:</p>
<pre>#look at coexistence between the two species within one cell
cell525&lt;-ldply(list_dat,function(x) c(x[525,&quot;N1&quot;],x[525,&quot;N2&quot;]))
cell525$Gen&lt;-1:61
cell525&lt;-melt(cell525,id.vars=&quot;Gen&quot;,
+ measure.vars=1:2,variable.name=&quot;Species&quot;,value.name=&quot;Pop&quot;)
ggplot(cell525,aes(x=Gen,y=Pop,color=Species))+geom_path()

#look at coexistence across the system
all&lt;-ldply(list_dat,function(x) c(sum(x[,&quot;N1&quot;]),sum(x[,&quot;N2&quot;])))
all$Gen&lt;-1:61
all&lt;-melt(all,id.vars=&quot;Gen&quot;,measure.vars=1:2,
+ variable.name=&quot;Species&quot;,value.name=&quot;Pop&quot;)
ggplot(all,aes(x=Gen,y=Pop,color=Species))+geom_path()</pre>
<p> </p>
<p><img alt="space1" class="wp-image-635" src="https://biologyforfun.files.wordpress.com/2016/06/space1.png?w=497"></p>
<p>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:</p>
<pre>#compare species-species abundance correlation at the beginning and at the end of the simulation
exmpl&lt;-rbind(list_dat[[1]],list_dat[[60]])
exmpl$Gen&lt;-rep(c(&quot;Gen : 1&quot;,&quot;Gen : 60&quot;),each=2500)
ggplot(exmpl,aes(x=N1,y=N2))+geom_point()+facet_grid(.~Gen)

<img alt="space3" class="wp-image-638" src="https://biologyforfun.files.wordpress.com/2016/06/space3.png?w=497">

#look at variation in spatial patterns at the beginning and the end of the simulation
exmplm&lt;-melt(exmpl,id.vars=c(&quot;x&quot;,&quot;y&quot;,&quot;Gen&quot;),
+ measure.vars=c(&quot;N1&quot;,&quot;N2&quot;),variable.name=&quot;Species&quot;,value.name=&quot;Pop&quot;)
ggplot(exmplm,aes(x,y,fill=Pop))+geom_raster()
 +scale_fill_viridis(option=&quot;plasma&quot;,direction=-1)+facet_grid(Gen~Species)</pre>
<p><img alt="space2" class="wp-image-640" src="https://biologyforfun.files.wordpress.com/2016/06/space2.png?w=497"></p>
<p> </p>
<p>Neat stuff, I am looking forward to implement other models in this book!</p>
<p> </p>
<p>Filed under: <a target="_blank" href="https://biologyforfun.wordpress.com/category/biological-stuff/">Biological Stuff</a>, <a target="_blank" href="https://biologyforfun.wordpress.com/category/r-and-stat/">R and Stat</a> Tagged: <a target="_blank" href="https://biologyforfun.wordpress.com/tag/ecology/">ecology</a>, <a target="_blank" href="https://biologyforfun.wordpress.com/tag/r/">R</a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/biologyforfun.wordpress.com/539/"><img alt="" src="http://feeds.wordpress.com/1.0/comments/biologyforfun.wordpress.com/539/" border="0"></a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/godelicious/biologyforfun.wordpress.com/539/"><img alt="" src="http://feeds.wordpress.com/1.0/delicious/biologyforfun.wordpress.com/539/" border="0"></a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gofacebook/biologyforfun.wordpress.com/539/"><img alt="" src="http://feeds.wordpress.com/1.0/facebook/biologyforfun.wordpress.com/539/" border="0"></a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gotwitter/biologyforfun.wordpress.com/539/"><img alt="" src="http://feeds.wordpress.com/1.0/twitter/biologyforfun.wordpress.com/539/" border="0"></a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gostumble/biologyforfun.wordpress.com/539/"><img alt="" src="http://feeds.wordpress.com/1.0/stumble/biologyforfun.wordpress.com/539/" border="0"></a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/godigg/biologyforfun.wordpress.com/539/"><img alt="" src="http://feeds.wordpress.com/1.0/digg/biologyforfun.wordpress.com/539/" border="0"></a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/goreddit/biologyforfun.wordpress.com/539/"><img alt="" src="http://feeds.wordpress.com/1.0/reddit/biologyforfun.wordpress.com/539/" border="0"></a> <img height="1" alt="" width="1" src="https://pixel.wp.com/b.gif?host=biologyforfun.wordpress.com&blog=42461605&%23038;post=539&%23038;subd=biologyforfun&%23038;ref=&%23038;feed=1" border="0"></p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="http://www.r-bloggers.com/exploring-spatial-patterns-and-coexistance/"></div></div><div><a data-url="http://www.r-bloggers.com/exploring-spatial-patterns-and-coexistance/" data-counturl="http://www.r-bloggers.com/exploring-spatial-patterns-and-coexistance/" data-count="vertical" data-text="Exploring Spatial Patterns and Coexistance" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://biologyforfun.wordpress.com/2016/06/11/exploring-spatial-patterns-and-coexistance/"> R – biologyforfun</a></strong>.</div>
<hr>
<a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="http://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="http://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="http://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="http://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="http://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="http://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="http://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="http://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="http://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="http://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="http://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="http://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="http://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="http://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="http://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="http://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/Atc7njoeFlA"></a></div>