Autómatas celulares (dinámica metapoblacional, sistema con competencia, sistema parásito-hospedero). Aplicación con R (paquete simecol)

Simulamos 3 casos: modelo metapoblacional (espacialmente implícito), modelo de competencia y modelo parásito-hospedero.

Modelo espacialmente implícito

#### Modelo espacialmente implícito: dinámica metapoblacional ##### producimos un dibujo con parches ocupados (color claro) y vacíos (color oscuro), de 10000 parches en un array de 100x100, pero el modelo no es espacialmente explícito, es solo una imagen. Las condiciones iniciales es de 100 parches ocupados aleatoriamente.
m=.15;e=.1;s=(1-e)
N<-matrix(rep(0,10000),nrow=100);xs<-sample(1:100);ys<-sample(1:100)
for (i in 1:100){N[xs[i],ys[i]]<-1 };image(1:100,1:100,N)
#queremos simulaciones de 1000 generaciones. Primero modelamos la supervivencia (o no) en parches ocupados. Cada celda en el universo toma un número aleatorio independiente desde una distribución uniforme (número real entre 0 y 1). Si el NA es <=s (=1-e) entonces el parche sobrevive otra generación. Si el NA es <s, el parche se extingue y N=0.
#calculamos la producción de propágulos (im) de los parches sobrevivientes. Asumimos que el asientamiento de los propágulos es aleatorio (caen en parches vacíos u ocupados).
for(t in 1:1000)
{
S <-matrix(runif(10000),nrow=100);N<-N*(S<s);
im<-floor(sum( N*m));placed<-matrix(sample(c(rep(1,im) ,rep (0,10000-im))),nrow=100)
N<-N+placed;N<-apply(N,2,function(x) ifelse(x>1,1,x)); image(1:100,1:100,N);box(col="red")
}
#como m>e la metapoblación persistirá. La solución analísitca a largo plazo es: 1-e/m=1-.1/.15/.333. En cualquier tiempo la proporción de ocupación actual es: sum(N)/length(N).


Modelo espacialmente explícito. Sistema de competencia (denso-dependencia local)

#### Modelo espacialmente explícito: DD local. competencia por sitios.
#Sean 2 especies que no pueden coexitir en un ambiente mezclado porque la fecundidad de la especie A es < que la fecundidad de la especie B, y esto conducirá tarde o temprano, en la exclusión competitiva de la especie B y la persistencia de la especie A. ej: insectos herbívoros o patógenos fúngicos.
#la sp B que se encuentra en la vecindad de la sp A impide el reclutamiento de la sp A cuando ésta se encuentra a más de un umbral T (individuos de sp A en el vecindario).
#idea: observar si la introducción de DD local en vecindario es suficiente para prevenir la exclusión competitiva y permitir la coexistencia a largo plazo de las 2 especies.
#bordes del universo:como todas las localizaciones neceitan tener el mismo número de vecinos en el modelo, utilizamos ‘wrap-around margins'
#definimos el vecindario, medinate una grilla cuadrada, una celda central y de 8 vecinos
plot(c(0,1),c(0,1),xaxt="n",yaxt="n",type="n",xlab="",ylab="");abline("v"=c(1/3,2/3));abline("h"=c(1/3,2/3))
xs<-c(.15,.5,.85,.15,.85,.15,.5,.85);ys<-c(.85,.85,.85,.5,.5,.15,.15,.15)
for (i in 1:8){ text(xs[i],ys[i],as.character(i)) }; text(.5,.5,"target cell")
#definimos los márgenes de las celdas. Sea el universo N de 100x100 celdas, la matriz de vecinos será de 102x102.
margins<-function(N)
{
edges<-matrix(rep(0,10404),nrow=102);edges[2:101,2:101]<-N;edges[1,2:101]<-N[100,];edges[102,2:101]<-N[1,];edges[2:101,1]<-N[,100]
edges[2:101,102]<-N[,1];edges[1,1]<-N[100,100];edges[102,102]<-N[1,1];edges[1,102]<-N[100,1];edges[102,1]<-N[1,100]
edges
}
#escribimos la funciónq ue cuenta el número de especies A en las 8 celdas vecinas de cada celda i,j:
nhood<-function(X,j,i) sum(X[(j-1):(j+1),(i-1):(i+1)]==1)
#seleccionamos los valores de los parámetros: tasa de reproducción de A y B, tasa de mortalidad de adultos y el número umbral T de especies A encima del cual no ocurre reclutamiento:
Ra<-3;Rb<-2.0;D<-0.25;s<-(1-D);T<-6
#condiciones iniciales: la mitad del universo con sp A y la otra mitad con sp. B
N<-matrix(c(rep(1,5000),rep(2,5000)),nrow=100);image(1:100,1:100,N)
#corremos la simulación 1000 pasos de tiempo
for (t in 1:1000)
{
S <-1*(matrix(runif(10000),nrow=100)<s) #para evaluar si un ocupante de una celda sobrevive o muere, comparamos un número aleatorio con distribución uniforme entre 0 y 1 con la tasa de sobrevivencia específica s=1-D. Si el NA<s el ocupante sobrevive, si NA>s muere.
N<-N*S;  space<-10000-sum(S)
nt<-margins(N) # calculamos la densidad de vecindario de A para cada celda
tots<-matrix(rep(0,10000),nrow=100)
for (a in 2:101) {
for (b in 2:101) {
tots[a-1,b-1]<-nhood(nt,a,b)
}}
seedsA<- sum(N==1)*Ra;  seedsB<- sum(N==2)*Rb #calculamos las semillas producidas por los sobrevivientes
all.seeds<-seedsA+seedsB;  fA=seedsA/all.seeds;  fB=1-fA
setA<-ceiling(10000*fA); placed<-matrix(sample(c(rep(1,setA) ,rep (2,10000-setA))),nrow=100) #las semillas se asientan aleatoriamente en el universo
#las semillas solo producen reclutas en celdas vacías (N[i,j]=0). Si el ganador de una celda vacía (placed) es la spB, toma la celda (if(placed[i,j]=2; N[i,j]<-2)). Si la spA gana la celda, necesitamos chequear que tiene <T vecinos de spA. Si lo hace la spA toma la celda. Si no lo hace la celda está dada a la spB (if(tots[i,j]>=T;N[i,j]<-2))
for (i in 1:100)
{
for(j in 1:100)
{
if(N[i,j] == 0 ){if(placed[i,j]== 2){N[i,j]<-2}} else {if(tots[i,j]>=T){N[i,j]<-2}else{N[i,j]<-1}}
}
}
image(1:100,1:100,N); box(col="red")  #dibujamos el mapa de spA en rojo y spB en blanco
}
#observamos que la spA auemnta en frecuencia a expensas de la sp.B. Eventualmente, sin emberago, la spA llega al pto donde su vecindario tiene muchos indiv de la sp.A, y su reclutamiento comienza a fallar. En el equilibrio, la spB persiste en celdas aisladas o en parches pequeños.
#si T=9 la sp.A conduce a la extinción de la sp.B
#si T=0 la sp.B toma el espacio


Modelo espacialmente explícito. Sistema parásito-Hospedero

#### Modelo espacialmente explícito: parásito-huésped.
#La interacción es inestable en modelos no-espaciales, con oscilaciones que aqumentan conduciendo a la rápida extinción de la sp.huesped y luego a la sp.parásito.
#en modelos espaciales (con movimiento a las 8 celdas vecinas), la interacción produce bonitos patrones espaciales que fluctúan con la abundancia de las 2 spp.
#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
}
##############################

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"
datos