Dinâmica de Metapopulações

<box 100% red|Antes de entrar em pânico…> Nesta prática iremos utilizar um programa que interage com você por meio de linhas de comando, o que significa que ele tem uma linguagem própria que vocês não devem conhecer. À primeira vista vai parecer assustador, mas vocês verão que não é tão ruim assim. Este roteiro atuará como seu guia e tradutor. Para interagir com o programa utilizaremos o GUI (graphic user interface), que apesar de básico, nesse caso ajuda muito. Sigam-no e ninguém sairá ferido. Bon voyage! </box>

Primeiro vamos abrir o tal programa e nos familiarizar um pouco com ele. Mas tente não perder muito tempo com essa primeira parte: se estiver demorando muito pule para a próxima seção.

Procure no desktop do seu computador um R azul e estiloso como o da figura acima. Se não acharem no desktop tentem em Iniciar → Todos os programas. Abram o programa e vocês verão uma tela árida e sem vida (Lembrem-se de que a Terra era assim também há alguns bilhões de anos atrás e veja no que ela se tornou). Note que há algumas poucas opções de comandos no menu superior.

Sem desanimar, vamos utilizar o R para fazer a coisa mais simples que ele sabe: operações aritméticas. Note que os comandos que devem ser inseridos no programa estarão sempre no formato de texto code, como abaixo. Não esqueça de que devem dar um ENTER depois de escrever cada linha de comando.

2+2
5-1
2*3
4/5

Notem que há um sinal de > em vermelho no canto esquerdo. Isso mostra o início da linha de comando. As linhas de comando, ou as ordens dadas ao programa, vão sempre ficar em vermelho, enquanto que as respostas que ele nos der vão sempre ficar em azul. Repare também que sempre, antes do resultado, aparece um [1]. Isso nada mais é do que um contador, quer ver? Vamos pedir agora uma seqüência de 1 a 100.

1:100

Note que o contador mostra a posição do primeiro valor de cada linha de resultado.

Continua... (opcional)

Nesta aula, iremos simular a dinâmica de uma metapopulação onde as taxas de extinção e de colonização são constantes (modelo de chuva de propágulos). Neste modelo há sempre uma fonte constante de imigrantes que podem colonizar qualquer mancha vazia. Relembrando, o modelo da dinâmica (variação da fração de manchas ocupadas no tempo) é descrito pela seguinte fórmula:

<m14> df/dt=C-E </m>

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

onde <m14>p_i</m> é a probabilidade de imigração ou colonização, <m14>p_e</m> é a probabilidade de extinção e <m14>f</m> é a fração de manchas ocupadas (número de manchas ocupadas / número total de manchas).

Em primeiro lugar, vamos estabelecer a probabilidade de colonização de manchas vazias (pi), a probabilidade de extinção em manchas ocupadas (pe) e a fração inicial de manchas ocupadas (fi) como 30%, 15% e 50%, respectivamente.

pi=0.3
pe=0.15
fi=0.5

Agora que já temos os parâmetros do nosso modelo (pi, pe e fi) vamos criar nossa paisagem virtual. Essa paisagem deve ser constituída por manchas de habitat que podem estar ocupadas ou não. Para tanto, criaremos uma matriz de 10 linhas e 10 colunas, sendo que cada célula dessa matriz representará uma mancha. Serão, portanto, 100 manchas.

paisag=array(0,dim=c(10,10,1))
paisag

Muito bem, mas temos somente a matriz no tempo inicial. Para acompanhar a dinâmica de ocupação de manchas é preciso criarmos uma terceira dimensão: o tempo. Por enquanto vamos criar 10 passos no tempo.

paisag=array(0,dim=c(10,10,10))
paisag

Ótimo! Agora temos nossa paisagem em 10 tempos diferentes.

Finalmente, antes de começarmos a brincadeira, precisamos definir quais manchas estarão ocupadas no tempo inicial. Vamos combinar que quando a mancha está ocupada ela recebe o valor 1 e quando está vazia recebe o valor 0. Do jeito que a matriz está, repleta de zeros, nenhuma mancha está ocupada. Então vamos preencher ao acaso algumas manchas, usando um fi de 50%:

paisag[,,1]=sample(c(rep(0,(100-fi*100)),rep(1,fi*100)))
paisag[,,1]

Agora sim, estamos prontos para a simulação.

O passo seguinte é fazer as coisas acontecerem! Mas para isso temos que ter em mente todas as possibilidades:

• Manchas ocupadas
  • podem permanecer ocupadas (1-pe = 0,85)
  • podem sofrer extinção local (pe = 0,15)
• Manchas vazias
  • podem permanecer vazias (1-pi = 0,7)
  • podem ser colonizadas (pi = 0,3)

Calma, não se assuste. O monstrinho abaixo vai tratar apenas das manchas que estavam ocupadas no tempo inicial.

paisag[,,2][paisag[,,1]==1]<-sample(c(0,1),sum(paisag[,,1]),replace=T,prob=c(pe,1-pe))
paisag[,,2]

Se você observar atentamente, algumas manchas que estavam ocupadas (1) ficaram vazias (0):

paisag[,,1] # paisagem no tempo inicial
paisag[,,2] # paisagem depois de um passo no tempo

Estes foram eventos de extinção local. Note que o número de manchas ocupadas diminuiu.

sum(paisag[,,1]) # total de manchas ocupadas no tempo inicial
sum(paisag[,,2]) # total de manchas ocupadas no tempo 2

Mas essa é apenas metade da história, não? Estamos esquecendo que as manchas vazias podem ser ocupadas. Vamos simular a colonização com um outro monstrinho:

paisag[,,2][paisag[,,1]==0]<-sample(c(0,1),10*10-sum(paisag[,,1]),replace=T,prob=c(1-pi,pi))
paisag[,,2]

Compare a paisagem no tempo inicial e no tempo 2 e veja se no final o número de manchas ocupadas aumentou ou diminuiu.

paisag[,,1] 
paisag[,,2]

sum(paisag[,,1]) # total de manchas ocupadas no tempo inicial
sum(paisag[,,2]) # total de manchas ocupadas no tempo 2

Note que algumas manchas que estavam ocupadas no tempo 1 ficaram vazias no tempo 2 e vice-versa. Note também que algumas manchas continuaram vazias e outras continuaram ocupadas. O próximo passo é calcular a fração de manchas ocupadas inicialmente (f1) e no tempo 2 (f2) e depois a diferença entre elas:

f1=sum(paisag[,,1])/100 # fração das manchas ocupadas no tempo 1
f2=sum(paisag[,,2])/100 # fração das manchas ocupadas no tempo 2

f1
f2
f2-f1

Eis a dinâmica de metapopulações! Mas até aqui vimos apenas a variação do f do tempo 1 para o tempo 2. Restam ainda 8 passos para chegarmos no tempo 10, só que para não termos que repetir todo o trabalho 8 vezes usaremos um comando chamado for, que fará o serviço por nós:

paisag=array(0,dim=c(10,10,10))
paisag[,,1]=sample(c(rep(0,(100-fi*100)),rep(1,fi*100)))
resultado=numeric()
for(t in 2:10){
	       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),10*10-sum(paisag[,,1]),replace=T,prob=c(1-pi,pi))
	       resultado[t-1]=sum(paisag[,,t])/(10*10)
	      }
resultado

Vôa lá! Agora temos as frações de manchas ocupadas (f) ao longo do tempo. Só para deixar nossos resultados mais bonitos, vamos colocá-los em uma tabela:

resultado=data.frame(t=1:10,f=c(fi,resultado))
resultado

Agora, para ficar mais bonito ainda, vamos criar um gráfico com esses resultados:

plot(1:10,resultado$f,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)

Como nossa filosofia de vida é (ou pelo menos deveria ser) melhorar sempre, vamos acrescentar a esse gráfico uma informação muito importante: a fração de manchas ocupadas no equilíbrio (F). O F é calculado da seguinte forma:

<m14>F=p_i/{p_i+p_e}</m>

Portanto…

F=pi/(pi+pe)
F

Agora que conhecemos o F vamos colocá-lo no gráfico na forma de uma linha horizontal:

plot(1:10,resultado$f,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)

O serviço sujo está feito, agora é hora da diversão! Como somos bonzinhos, resolvemos criar uma função para facilitar a vida de vocês. Com essa função vocês poderão variar os parâmetros do nosso modelo à vontade, sem medo de ser feliz. Por favor, retribuam a gentileza e testem vários valores para cada parâmetro e vejam o que acontece com nossas metapopulações virtuais. Abaixo a função:

metapop=function(tf,c,l,fi,pi,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){
	       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=pi/(pi+pe)

	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)
	}

Maravilha! A função acima tem 6 parâmetros (tf, c, l, fi, pi e pe), mas para que a função funcione é preciso atribuir um valor para cada parâmetro, como no exemplo abaixo:

metapop(tf=100,c=10,l=10,fi=0.5,pi=0.3,pe=0.15)

<box 50%|Parâmetros da função metapop>

• tf - número total de passos no tempo
• c - número de colunas da matriz paisag
• l - número de linhas da matriz paisag
• fi - fração de manchas ocupadas inicial
• pi - probabilidade de colonização
• pe - probabilidade de extinção

</box> Agora é com vocês! Tentem responder às seguintes perguntas:

• O que acontece quando o número de manchas é muito grande? E quando é muito pequeno?
• E quando a fração inicial (fi) é muito diferente da fração em equilíbrio (F)?
• Como a população pode se extinguir da paisagem?
• O que acontece quando as probabilidades de extinção e de colonização são muito altas?
• E quando elas são muito baixas?
• O que altera a frequência da trajetória de f? E a amplitude?
• Do que depende a velocidade com que a trajetória se aproxima do F?
• Em que situações a trajetória chega a coincidir exatamente com a linha do F?
• Quando você rodar várias mais de uma vez a função com os mesmos parâmetros o resultado é o mesmo? Por quê?

Para vocês que chegaram vivos até aqui, uma recompensa. Rode os comandos abaixo e descubra:

metapop2=function(tf,c,l,fi,pi,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){
	       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)
	      }
	return(paisag)
	}


anima=function(tf,c,l,fi,pi,pe){
	dados=metapop2(tf,c,l,fi,pi,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)
	}
	}

anima(tf=25,c=10,l=10,fi=.1,pi=0.1,pe=0.1)

Aqui você também pode mexer nos parâmetros e ver, literalmente, o que acontece:

anima(tf=25,c=100,l=100,fi=.01,pi=0.2,pe=0.5)

E aí, gostaram?

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)

Script das funções:

• metapop.pdf