Podemos eliminar do modelo anterior o pressuposto de uma chuva de propágulos constante e fazer com que a colonização seja uma função do número de lugares ocupados. Em uma formulação simples desse modelo, a fonte de propágulos é unicamente interna (sistema fechado) e a probabilidade de colonização varia de forma linear à proporção de lugares ocupados.

Dessa forma, nosso modelo não terá mais uma probabilidade de colonização constante (pi), mas sim uma probabilidade de colonização dependente do número de manchas ocupadas:

<m14>p_i=if</m> ; onde i é uma constante que indica quanto aumenta a pi a cada nova mancha que é ocupada.

Portanto, quanto mais manchas ocupadas, maior a chance de colonização das manchas vazias. Substituindo pi na equação antiga temos:

<m14>df/dt=if(1-f)-p_e f</m>

O cálculo da fração de manchas ocupadas no equilíbrio (F) também é modificado para:

<m14>F=1-p_e /i</m>

meta.inter=function(tf,c,l,fi,i,pe){
	paisag=array(0,dim=c(l,c,tf))
	paisag[,,1]=sample(c(rep(0,round(c*l-fi*c*l)),rep(1,round(fi*c*l))))
	resultado=numeric()
	for(t in 2:tf){
		 pi=i*sum(paisag[,,t-1])/(c*l)
	       paisag[,,t][paisag[,,(t-1)]==1]<-sample(c(0,1),sum(paisag[,,1]),replace=T,prob=c(pe,1-pe))
	       paisag[,,t][paisag[,,(t-1)]==0]<-sample(c(0,1),c*l-sum(paisag[,,1]),replace=T,prob=c(1-pi,pi))
	       resultado[t-1]=sum(paisag[,,t])/(c*l)
	      }

	F=1-(pe/i)

	plot(1:tf,c(fi,resultado),type="l",xlab="Tempo",ylab="Fração de manchas ocupadas",
	ylim=c(0,1),main="Dinâmica de ocupação de manchas",font.lab=2,lwd=2)
	abline(h=F,col=2,lwd=2,lty=2)
	
      return(paisag)
	}

meta.inter(tf=100,c=10,l=10,fi=.1,i=1,pe=0.5)

Brinque um pouco com o modelo fazendo variar os parâmetros do modelo e pense nas seguintes perguntas:

• A posição de uma mancha na paisagem influencia a pi e a pe dessa mancha? Qual seria um modelo mais realista?
• Neste novo modelo as extinções regionais são mais comuns? Por quê?
• Por que há certas combinações de i e pe que não podem existir?
• Qual o significado de um F negativo?
• Você consegue perceber alguma diferença nos resultados dos dois modelos?

Para finalizar, uma última animaçãozinha:

meta.inter2=function(tf,c,l,fi,i,pe){
	paisag=array(0,dim=c(l,c,tf))
	paisag[,,1]=sample(c(rep(0,round(c*l-fi*c*l)),rep(1,round(fi*c*l))))
	resultado=numeric()
	for(t in 2:tf){
		 pi=i*sum(paisag[,,t-1])/(c*l)
	       paisag[,,t][paisag[,,(t-1)]==1]<-sample(c(0,1),sum(paisag[,,1]),replace=T,prob=c(pe,1-pe))
	       paisag[,,t][paisag[,,(t-1)]==0]<-sample(c(0,1),c*l-sum(paisag[,,1]),replace=T,prob=c(1-pi,pi))
	       resultado[t-1]=sum(paisag[,,t])/(c*l)
	      }

	F=1-(pe/i)

      return(paisag)
	}

anima2=function(tf,c,l,fi,i,pe){
	dados=meta.inter2(tf,c,l,fi,i,pe)
	for(i in 1:tf){
	image(dados[,,i], main=("Ocupação de manchas"),col=c("white","red"),bty="n",xaxt='n',yaxt='n')
	grid(c,l)
	Sys.sleep(.5)
	}
	}

anima2(tf=25,c=10,l=10,fi=.1,i=1,pe=0.1)

Abaixo o script das funções. A extensão está como pdf, mas na realidade é um arquivo de texto simples. Para abrir deve salvar em seu computador e modificar a extensão para .txt. Depois disso pode abrir em qualquer editor de texto (Word, Notepad, ou com o próprio editor do R)

• metapop.pdf