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 caso há sempre uma fonte constante de imigrantes que podem colonizar qualquer mancha vazia.

Primeiro vamos conhecer um novo comando do R que cria uma matriz de dados tridimensional. O comando é array. Teste a seguinte linha para ver um array de 3 dimensões.

array(NA,dim=c(4,4,3))

Veja que o R criou 3 matrizes de 4 linhas por 4 colunas. Essa matriz assemelha-se a um bolo com 3 camadas, porém quando o objeto é chamado o R mostra as camadas separadamente. No exercício abaixo a primeira dessas 3 matrizes representará o tempo t = 1, e cada célula da matriz representará uma mancha no espaço, que poderá ou não estar ocupada.

• Crie uma função chamada metapop com os seguintes argumentos:
  • nl: número de linhas;
  • nc: número de colunas;
  • tmax: número máximo de vezes que você deseja simular a dinâmica de metapopulações;
  • psi0: proporção de manchas ocupadas no início da simulação;
  • cln: probabilidade de colonização de manchas vazias;
  • ext: probabilidade de extinção local em manchas ocupadas
• O número de linhas nl e o número de colunas nc determinarão quantas manchas potencialmente ocupáveis haverá na metapopulação (nc*nl = número de manchas). Use os seguintes comandos para começar a sua função:

  ## crie esses objetos para ir testando as linhas da sua função enquanto cria
tmax = 20; nc = 10; nl = 10; psi0 = 0.5; cln = 0.3; ext = 0.15  
  ## Para criar a função.
metapop<-function(nl,nc,tmax,psi0,cln,ext){   
  ## Isso vai criar tmax matrizes com nl linhas e nc colunas.
resulta<-array(0,dim=c(nl,nc,tmax))  
  ## como vamos usar muito o nl*nc (número de manchas), vamos criá-lo

Agora, para cada mancha você deverá jogar uma moeda para decidir se a mancha está ou não ocupada no tempo t =0. A probabilidade de estar ocupada deve ser psi0 (a probabilidade de ocorrência inicial). Para isso use a função sample().
Usando sample(c(0,1),1) obtemos uma amostra de um valor que pode ser 0 ou 1, com igual probabilidade. Mas nós precisamos amostrar nl*nc valores de 0 ou 1, pois temos esse número de manchas.

Então vamos usar:

sample(c(0,1),nl*nc,replace=T )

onde replace=T significa amostragem com reposição.

Porém, ainda falta informarmos a probabilidade de ocorrência. Se amostrarmos com igual probabilidade de obter 0 (mancha desocupada) ou 1 (mancha ocupada) a probabilidade de ocorrência será sempre 0.5, mas não é isso que pretendemos. Queremos ter a opção de alterar a probabilidade de ocorrência. Para isso usamos o argumento “prob=c(1-psi0,psi0)” que especifica a probabilidade de obter os valores indicados dentro da função sample. Ou seja, 1-psi0 será a probabilidade de um 0 ser amostrado e psi0 será a probabilidade de um 1 ser amostrado.

Então vamos simular a ocupação no tempo t=0 de acordo com psi0:

resulta[ , ,1]<-sample(c(0,1),nl*nc,replace=T,prob=c(1-psi0,psi0))

A probabilidade de ocorrência psi0 irá variar nos próximos passos, então crie um objeto chamado de resultados “res.psi” e salve o psi0 na posição 1 (psi0 é o valor inicial de fs). No Gotelli f é a fração de manchas ocupadas (número de manchas ocupadas / número de manchas).

Use um for() para simular o que acontece com a metapopulação em cada tempo subseqüente, até que tmax seja atingido. Repare que já preencheu o array de resultado na primeira posição relativa ao tempo inicial

Agora precisamos criar uma condição para saber se uma mancha que estava ocupada continuará ocupada e se uma mancha que estava desocupada continuará desocupada, ou se será ocupada.

Veja o que apareçe no objeto resulta[,,1]: apenas valores 0 e 1, significando ausência ou ocupação de mancha no tempo t=1. No tempo t+1 ou no caso t=2 as possibilidades são:

• Manchas desocupadas (valores 0) podem ser colonizados e passar a 1, ou continuar 0.
• Manchas ocupadas (valores 1)podem ter extinção e passsar a 0, ou continuar 1

Para o primeiro caso a probabilidade de 0 passar a 1 é cln e a probabilidade de 1 passar a 0 é ext. Uma forma de fazer isso é:

  num1=sum(resulta[,,t-1])
  num0=(nl*nc)-num1
  resulta[,,2][resulta[,,1]==1]<-sample(c(0,1),num1,replace=T,prob=c(ext,1-ext))
  ...

Esse comando pode ser lido da seguinte forma: nas posições do resulta t2 onde resulta t1 são iguais a 1, coloque 0 ou 1 com probabilidade ext e 1-ext respectivamente.

Faça o mesmo para as posições em t1 iguais a 0, pensando na colonização. Pronto, agora e só mudar o valores de tempo pelo seu indexador para que em cada ciclo um novo tempo seja calculado. Para finalizar o ciclo só falta calcular o res.psi de cada tempo.

res.psi[t]=sum(resulta[,,t])/(nl*nc)

<box 100% |Dica avançada> Sugerimos que use a solução anterior por ser mais intuitiva. Caso queira se aprofundar nos comandos do R, veja a solução abaixo para o mesmo problema: substituir os dados de ocupação de manchas de um tempo para outro, utilizando as probabilidades de extinção e colonização:

É possível fazer o mesmo com uma única linha de comando utilizando a função ifelse(). Ela funciona assim:

ifelse(condição, “faça isso”, “se não faça aquilo”)

Então tente algo assim:

 resulta[,,t]<-ifelse(Passo 1, Passo 2, Passo3) 

Veja abaixo as instruções para os Passos 1, 2 e 3.
Instruções para os passos 1, 2 e 3 acima.
Na função ifelse acima você deverá, no Passo 1, criar uma condição para saber se uma determinada mancha estava ou não ocupada no tempo t-1.
Caso a mancha estivesse desocupada você deverá, no Passo 2, dizer o que acontecerá com ela no tempo t. Para isso gere um número aleatório de uma distribuição uniforme e veja se esse valor é menor que a taxa de colonização, se o valor for menor a mancha deverá ser ocupada, se for maior ela permanece desocupada. Vamos ver como fazer isso:
Atenção, as duas linhas abaixo são para você testar no R e entender os comandos, estas linhas não devem estar no seu script final.
runif(1)<cln Veja que o resultado será TRUE ou FALSE, mas nós queremos que seja 0 ou 1

Precisamos transformar o TRUE e o FALSE em valores numéricos, usando a função as.numeric()

as.numeric(runif(1)<cln) 	

veja que o resultado agora será 0 ou 1.
Agora você já sabe como fazer o Passo 2, mas lembre que no runif(1'') você não deverá ter apenas 1 valor, e sim nc*nl valores, ou seja, um valor pra cada mancha.

Passo 3: Aprenda com o Passo 2.

Fim das instruções para os passos 1, 2 e 3.
</box>

Agora faça um gráfico, dentro da sua função, que mostre a trajetória da proporção de manchas ocupadas ao longo do tempo. Acrescente uma reta mostrando o valor da fração de manchas ocupadas no equilíbrio f.equilibrio.
Equação 4.4 do Gotelli: f=Pi/(Pi+Pe)

Primeiro calcule o f.equilibrio

f.equilibrio<- ??

plot(1:tmax),fs ,type="l",xlab="tempo",ylab="fs",ylim=c(0,1),main="Metapopulação ")
abline(h=f.equilibrio, col=2)

Agora vamos criar uma animação para olhar nossos resultados ao longo do tempo. Este mostra uma animação da dinâmica da metapopulação ao longo do tempo.
Primeiro salve o resultado de uma simulação da função metapop e crie um gráfico do tempo 1.

   image(dados[,,i], main=("Ocupação de manchas"),col=c("white","red"))
   grid(nl,nc)
 

Agora é só criar o ciclo com for para que os gráficos de cada tempo sejam sobrepostos no mesmo dispositivo gráfico. Para isso vamos criar uma função anima.meta:

1. O único argumento é o dado que será igual ao resultado da função metapop.
2. crie um ciclo de 1 até tmax que é igual ao tamanho da dimensão três do seus dados
3. faça com que em cada ciclo seja feito um gráfico usando image dos dados do tempo t
4. crie um lapso de tempo para que os gráficos possam ser visualizados antes de trocar, colocando o comando abaixo:

Sys.sleep(.5)

• pronto!

Use sua função e altere os parâmetros iniciais para ver se a metapopulação persiste, dependendo das condições iniciais. Se tiver tempo, procure alterar a sua função para implementar os outros modelos descritos na aula téorica.

<box 90% blue| Uma animação mais informativa> Uma forma de representar a dinâmica espacial da metapopulação é incluir no mapa não só os espaços ocupados e vagos mas entre esses aqueles que foram recentemente colonizados ou extintos. Ou seja, teremos quatro tipos de informações em cada momento: ocupado antigo, ocupado recente, vago, extinção recente.

meta.anima2=function(dados)
{
nsim=dim(dados)[3]
image(dados[,,1], main=("Ocupacao de manchas"),col=c("white","green"),sub="tempo=1")
	for(i in 2:nsim)
	{
	conta12=dados[,,(i-1)]+ (2*dados[,,i])
	image(conta12, main=("Ocupacao de manchas"),col=c("white", "red","green",
"darkgreen"),sub=paste("tempo= ",i))
	Sys.sleep(.5)
	}
}

</box>

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. Nossa função começa então diferente com relação aos argumentos necessários:

• sai o argumento cln;
• entra o argumento ef: quanto a proporção de áreas ocupadas afeta a colonização.

meta.inter<-function(nl,nc,tmax,psi0,ext, ef)
{   
...

• calcule a probabilidade de colonização no inicio de cada ciclo:

...
  for(t in 2: tmax)
  {
   cln=ef*res.psi[t-1]
   ...
   ...
  } 

• modifique o valor da proporção de áreas ocupadas no equilíbrio para:

...
f.equilibrio=1-(ext/ef)
...   

• com essas mudanças a novo modelo deve funcionar!