library(spdep) #Required to use the function cell2nb() N<-25 #Number of rows and columns in the grid (not the same as N in the paper) time.steps<-2500 #Number of time steps amax<-60 #Maximum time steps F. cupressoides can live bmax<-30 #Maximum time steps Nothofagus spp. can live disturbed.ratio<-0.5 #Probability a cell is initially disturbed F.mature<-7 #Age of maturity for F. cupressoides N.mature<-2 #Age of maturity for Nothofagus spp. catastrophe.prob<-0.1 #Probability a catastrophe happens in a cell at each time step neighbor.indices<-cell2nb(N,N,type="rook",torus=TRUE) #creates a list such that each element #stores the indices of the neighbors of that element F.pop<-as.numeric(runif(N*N)<0.3) #Is the cell occupied by F. cupressoides? N.pop<-as.numeric(runif(N*N)<0.5) #Is the cell occupied by Nothofagus spp.? s<-as.numeric(runif(N*N)0))/(N*N) #Initial density of F. cupressoides meanN.density[1]<-sum(as.numeric(b>0))/(N*N) #Initial density of Nothofagus spp. for(i in 1:time.steps){ catastrophe<-as.numeric(runif(length(s))F.mature))/4)) == 1){ a[j]<-1 timea.colonized[j]<-1 } else{ } } else if(a[j] < amax){ #if cell is colonized by F. cupressoides then add to occupation time a[j]<-a[j]+1 timea.colonized[j]<-timea.colonized[j]+1 } else{ #if occupation time has reach amax then F. cupressoides dies a[j]<-0 timea.colonized[j]<-0 } } else{ #If soil is not disturbed... if(b[j] == 0){ #If cell is not colonized by Nothofagus spp. then chck if it becomes colonized if(as.numeric(runif(1)<(sum(as.numeric((b[neighbor.indices[[j]]])>F.mature))/4)) == 1){ b[j]<-1 } else {} } else if(b[j] < bmax){ #If cell is colonized by Nothofagus spp. then add to occupation time b[j]<-b[j]+1 } else{ #If occupation time has reach bmax then keep it at bmax a[j]<-bmax } if(a[j] == 0 && b[j] == 0){ #If cell is not colonized by either species check if F. cupressoides colonizes if(as.numeric(runif(1)<(sum(as.numeric((a[neighbor.indices[[j]]])>F.mature))/4)) == 1){ a[j]<-1 timea.colonized[j]<-1 } else{} } #If cell is colonized by both species then add to occupation time of F. cupressoides else if(a[j] < amax && b[j] < bmax){ a[j] < a[j] + 1 timea.colonized[j]<-timea.colonized[j]+1 } else if(a[j] == amax || b[j] == bmax){ #If either species has reach its max occupation time then F. cupressoides dies a[j]<-0 timea.colonized[j]<-0 } else{ } } if(a[j] == 15){ #If F. cupressoides occupation time is 15 then soil in cell is no longer disturbed s[j]<-1 } else { } } } meanF.density[i+1]<-sum(as.numeric(a>0))/(N*N) #Store density of F. cupressoides meanN.density[i+1]<-sum(as.numeric(b>0))/(N*N) #Sotre density of Nothofagus spp. } plot(c(0,time.steps),c(0,1),type="n",main="Population Density", xlab="time",ylab="Density") points(0:time.steps,meanF.density) points(0:time.steps,meanN.density,pch=2) legend(1700,0.3,legend=c("F. cupressoides", "Nothofagus spp."), pch=c(1,2))