# A script that explores some properties of
# a metapopulation in which colonization
# only comes from within the metapopulation
# model with parameters explained below

pops=10	 # Number of population units
e=.1; c=.2;	 # extinction & colonization rates
peq=1-e/c	 # Equilibrium value for the
		 #   proportion of patches occupied
yrs=100	 # Years to simulate

Ns=matrix(NaN,yrs,pops)    # Matrix with 0/1 for each
Ns[1,]=rbinom(pops,1,peq)  #  population in each year

for (i in 1:(yrs-1)) {
   for (j in 1:pops)
	if (Ns[i,j]==0) Ns[i+1,j]=rbinom(1,1,c*sum(Ns[i,]/pops))
      else Ns[i+1,j]=rbinom(1,1,(1-e))
}

# Proportion of all populations occupied in each year
Frac=apply(Ns,1,sum)/pops   

plot(1:yrs,Frac,ylim=c(0,1),'o')
abline(h=peq,col='green')
mtext(paste("Avg Ppn Occupied = ", mean(Frac), sep = ""),
side = 3, line = 1, cex = 1, font = 2)
mtext(paste("Expected Ppn Occupied = ", mean(peq), sep = ""),
side = 3, line = 3, cex = 1, font = 2)

# print key features to screen
mean(Frac)
peq