Automatas celulares con R.

Les dejo algunos ejemplos de autómatas celulares aplicados en dinámica poblacional, ya iré actualizando la página con más scripts.
También pueden mirar el excelente trabajo de Ginger Booth: http://www.gingerbooth.com/pages/demosecology.html y http://desert.bgu.ac.il/desert/EngSite.aspx?SiteId=3354&ItemId=3365.

Ejemplo 1. Modelo espacial de Nicholson-Bailey.
Construimos una grilla espacial de interacción entre parásitos y hospederos y modelamos su dinámica poblacional.
Incluiré próximamente el contenido teórico de este tema, pero mientras si quieren puede mirar los siguientes archivos: simecol y cómo hacerlo y más sobre el paquete

################################################
#Modelo Nicholson-Bailey: función host para calcular la próxima población de hospedadores como función de el número actual de hospedadores y parásitos (N y P), y otra función parasite para calcular la p´roxima población de parásitos como función de N y P
################################################
 ##############################

#función de Nicholson-Bailey
host=function(N,P,r,a) N*exp(r-a*P)
parasite=function(N,P,a) N*(1-exp(-a*P))

#definiciones de los bordes
host.edges=function(N,Nrow)
 {
 Hedges=matrix(rep(0,(Nrow+2)^2),nrow=(Nrow+2))
 Hedges[2:(Nrow+1),2:(Nrow+1)]=N
 Hedges[1,2:(Nrow+1)]=N[Nrow,]
 Hedges[(Nrow+2),2:(Nrow+1)]=N[1,]
 Hedges[2:(Nrow+1),1]=N[,Nrow]
 Hedges[2:(Nrow+1),(Nrow+2)]=N[,1]
 Hedges[1,1]=N[Nrow,Nrow]
 Hedges[(Nrow+2),(Nrow+2)]=N[1,1]
 Hedges[1,(Nrow+2)]=N[Nrow,1]
 Hedges[(Nrow+2),1]=N[1,Nrow]
 Hedges
 }

parasite.edges=function(P,Nrow)
 {
 Pedges=matrix(rep(0,(Nrow+2)^2),nrow=(Nrow+2))
 Pedges[2:(Nrow+1),2:(Nrow+1)]=P
 Pedges[1,2:(Nrow+1)]=P[Nrow,]
 Pedges[(Nrow+2),2:(Nrow+1)]=P[1,]
 Pedges[2:(Nrow+1),1]=P[,Nrow]
 Pedges[2:(Nrow+1),(Nrow+2)]=P[,1]
 Pedges[1,1]=P[Nrow,Nrow]
 Pedges[(Nrow+2),(Nrow+2)]=P[1,1]
 Pedges[1,(Nrow+2)]=P[Nrow,1]
 Pedges[(Nrow+2),1]=P[1,Nrow]
 Pedges
 }
#función de vecindario
nhood=function(X,j,i) sum(X[(j-1):(j+1),(i-1):(i+1)])

#función de llegada de migrantes a cada celda
h.migration<-function(Hedges,Nrow) #The number of host migrants arriving in every cell is calculated as follows:
 {
 Hmigs<-matrix(rep(0,Nrow^2),nrow=Nrow)
 for (a in 2:(Nrow+1))
   {
 for (b in 2:(Nrow+1))
     {
     Hmigs[a-1,b-1]<-nhood(Hedges,a,b)
     }
   }
 Hmigs
 }

 p.migration<-function(Pedges,Nrow)  # The number of parasites migrants is given by:
 {
 Pmigs<-matrix(rep(0,Nrow^2),nrow=Nrow)
 for (a in 2:(Nrow+1))
   {
   for (b in 2:(Nrow+1))
     {
     Pmigs[a-1,b-1]<-nhood(Pedges,a,b)
     }
   }
 Pmigs
 }
 ##############################

#### modificado #################################
parasito_hospedero<-function(Tiempo=100,Xinitial=33,Yinitial=33,Nrow=100,NN=2,NP=2,r=.4,a=.1,Hmr=.1,Pmr=.9)
{

  ##seleccionamos los valores de los parámetros: dinámicas del hospedador (r) y parásito (r), tasas de migración Hmr y Pmr: en este caso los hospedadores son relativamente sedentarios y el parásito es altamente móvil
N=matrix(rep(0,Nrow^2),nrow=Nrow) #matrices de abundancias del hospedador (N) y el parásito (P)
P=matrix(rep(0,Nrow^2),nrow=Nrow)
N[Xinitial,Yinitial]=NN #condiciones iniciales 200 hospedadores y Nrow parásitos en una única calda de localización (33,33)
P[Xinitial,Yinitial]=NP
Nt<-NN;Pt<-NP
par(mfrow=c(2,1));  image(N+P) #eliminar si deseo obtener solo el gráfico de Nt y Pt vs tiempo
plot(NN~1, xlim=c(0,Tiempo),ylim=c(0,NN*2),pch=20,col="red"); points(NP~1,pch=20,col="blue")
 for (t in 2:Tiempo) #The simulation begins here, and runs for 600 generations:
 {
 he<-host.edges(N,Nrow=Nrow)
 pe<-parasite.edges(P,Nrow=Nrow)
 Hmigs<-h.migration(he,Nrow=Nrow)
 Pmigs<-p.migration(pe,Nrow=Nrow)
 N<-N-Hmr*N+Hmr*Hmigs/9
 P<-P-Pmr*P+Pmr*Pmigs/9
 Ni<-host(N,P,r,a)
 P<-parasite(N,P,a)
 N<-Ni
 Nt<-c(Nt,sum(N))
 Pt<-c(Pt,sum(P))
 image(N+P);    TT<-1:t;plot((Nt/100)~TT,xlim=c(0,Tiempo),ylim=c(0,NN*2),pch=20,type="o",col="red");points((Pt/100)~TT,pch=20,type="o",col="blue") #eliminar si deseo obtener solo el gráfico de Nt y Pt vs tiempo
 #image(1:Nrow,1:Nrow,N)
#TT<-1:t;points((Nt/100)~TT,pch=20,type="o",col="red");points((Pt/100)~TT,pch=20,type="o",col="blue") #agregar si deseo obtener el gráfico de Nt y Pt vs Tiempo
 }
out<-list(Hospederos=Nt,Parasitos=Pt)
out
}

datos<-parasito_hospedero(Tiempo=100,Xinitial=33,Yinitial=33,Nrow=100,NN=200, NP=100, r=.4,a=.1,Hmr=.1,Pmr=.9) #es el ejemplo del libro "R book"