A Teoria Neutra é um modelo de processos estocásticos de nascimentos, mortes, especiações e migrações. As probabilidades de cada um destes eventos ocorrerem definem uma dinâmica surpreendente. A melhor maneira de entender isto é simular este processo, como faremos nos exercícios a seguir.

Um processo estocástico acontece quando temos mais de um estado possível para um sistema, e podemos pular para cada um com uma certa probabilidade. A Teoria Neutra se encaixa em uma classe particular destes processos, chamada Cadeia de Markov. Nessa Cadeia de Markov o tempo é discreto, e a cada intervalo o sistema pode mudar de estado, com uma certa probabilidade. As probabilidades de mudança de um estado para outro dependem apenas do estado presente ¹⁾.

Para entender duas propriedades importantes destes processos, vamos simular dois casos simples.

Imagine um bêbado andando sempre para frente em uma enorme planície, mas que tem um abismo em um dos lados. A cada passo para frente, ele cambaleia um certo numero de passos para a direção do abismo ou da planície, com igual probabilidade.

Este é um dos processos Markovianos mais simples, chamado caminhada aleatória (random walk) em uma dimensão ²⁾. Se o bêbado cai no abismo a caminhada acaba (e o bêbado também), uma condição que chamamos de fronteira de absorção (absorbing boundary).

O que podemos prever deste processo? Vamos soltar alguns bêbados neste mundo virtual. Para isto usaremos uma função em R, cujo código segue abaixo. Copie e cole este código na linha de comando do R.

bebado <- function(n=1,step=1,ciclo=1e5,cont=1e3,x1=NULL){
  if(is.null(x1)){
    x1 <- sample(1:200,n,replace=TRUE)
  }
  results <- matrix(NA,nrow=1+ciclo/cont,ncol=n) 
  results[1,] <- x1
  X <- x1
  for(i in 2:(1+ciclo/cont)){
    for(j in 1:cont){
      X[X<=0] <- NA
      X <- X +sample(c(step,-1*step),n,replace=T)
    }
    results[i,] <- X
  }
  results[is.na(results)] <- 0
  time <- seq(0,ciclo,by=cont)
  matplot(time,results,type="l", col=rainbow(n),lwd=2, xlab="Passos", ylab="Distância do Abismo")
  abline(h=0,lwd=4)
}

Esta função tem quatro argumentos:

• n: Número de bêbados
• step: número de passos para o lado que cada bêbado dá a cada instante de tempo
• ciclo: tempo total da simulação (no caso é o mesmo que o número de passos para frente).
• cont: intervalo de registro dos dados. Este argumento serve para poupar memória, já que guardar todos os passos em simulações longas pode deixar a simulação muito lenta. Mantendo o valor default a posição de cada bêbado é registrada a cada 1000 intervalos de tempo.

A função sorteia a posição inicial dos bêbados entre 1 e 200 passos de distância do abismo. Vamos soltar dez bêbados, que cambaleiam 10 passos a cada intervalo, por dez mil intervalos de tempo:

bebado(n=10,step=10,ciclo=1e4,cont=1e3)

Repita esta simulação algumas vezes, e tente identificar um padrão.

O que acontece se deixamos os bêbados um pouco menos cambaleantes? Experimente reduzir para dois os passos laterais:

 
set.seed(42) # semente de números aleatórios
bebado(n=10,step=2,ciclo=1e4,cont=1e3)

Estes bêbados que balançam menos estão menos sujeitos a terminar no abismo? Certifique-se disto aumentando o número de intervalos de tempo:

set.seed(42)
bebado(n=10,step=2,ciclo=5e4,cont=1e3)

Esta caminhada aleatória com fronteira de absorção tem um único desfecho, dado tempo suficiente. Qual é?

Agora vamos imaginar um jogo de apostas entre dois jogadores, sem empates. A cada rodada o perdedor da aposta paga uma quantia fixa ao ganhador. Os dois jogadores tem a mesma probabilidade de ganhar a cada rodada. Esta é uma dinâmica de soma zero, pois o valor total em jogo não se altera. O que muda é apenas a fração deste total em poder de cada jogador.

Em nossa simulação, o jogo só termina quando acaba, ou seja, quando um dos dois jogadores perde todo o dinheiro. Vamos simular esta situação com a função em R abaixo:

jogo <- function(aposta=1,total=20){
  X <- total/2
  results <- X
  while(X>0&X<total){
    X <- X+sample(c(aposta,-1*aposta),1)
    results <- c(results,X)
  }
  plot(1:length(results),results, type="l", col="blue",ylim=c(0,total), xlab="N de rodadas", ylab="Valor")
  lines(1:length(results),total-results, col="red")
  abline(h=c(0,total),lty=2)
}

Esta função tem apenas dois argumentos:

• Aposta: que é o valor pago pelo perdedor a cada aposta
• Total: o valor total de dinheiro em jogo. No início este total é dividido igualmente entre os jogadores.

A simulação transcorre até o final do jogo. Investigue o efeito do aumento do total em jogo sobre o tempo para que o jogo acabe.

Este jogo é também um processo de caminhada aleatória em uma dimensão. Explique porque.

Agora que entendemos algumas propriedades básicas de cadeias Markovianas simples vamos construir o modelo estocástico da Teoria Neutra, passo a passo.

Vamos começar com um modelo para a comunidade em um dado local, usando um jogo de soma zero, similar ao jogo de apostas da seção anterior. As regras são:

1. A comunidade tem um total fixo de indivíduos, $$J$$
2. A cada intervalo de tempo, um dos indivíduos é sorteado para morrer
3. Em seguida, os indivíduos remanescentes são sorteados, para definir quem produzirá o filhote que ocupará o lugar do morto.

Para simular este processo, temos mais uma super-função em R logo abaixo. Copie-a e cole-a na linha de comando para usá-la:

rich <- function(x)length(unique(x)) ## funcao auxiliar

sim.hub1=function(S= 100, j=10, D=1, ciclo=2e4, step=1000){ 
  ## Tamanho da comunidade
  J <- S*j
  ##Matrizes para guardar os resultados
  ## matriz da especie de cada individuo por ciclo
  ind.mat=matrix(nrow=J,ncol=1+ciclo/step) 
  ##CONDICOES INICIAIS##
  ##Deduzidas de acordo com o modelo de Hubbell:
  ## Todas as especies comecam com o mesmo numero de individuos (j=J/S)
  ind.mat[,1] <- rep(1:S,each=j)
  cod.sp <- ind.mat[,1]
  ##Aqui comecam as simulacoes
  for(i in 2:(1+ciclo/step)){
    for(j in 1:step){
      ##Indice dos individuos que morrem
      morte <- sample(1:J,D)
      ##Indice dos individuos que produzem filhotes para substituir os mortos
      novos <- sample(1:J,D,replace=T)
      ##Substituindo
      cod.sp[morte]<-cod.sp[novos]
    }
    ## A cada step ciclos os resultados sao gravados
    ind.mat[,i] <- cod.sp
  }
  tempo <- seq(0,ciclo,by=step)
  colnames(ind.mat) <- tempo
  invisible(ind.mat)
  plot(tempo,apply(ind.mat,2,rich), xlab="Tempo (ciclos)", ylab="N de espécies", type="l",
       main=paste("Dinâmica Neutra sem Colonização", "\n S=",S," J=",J)
       ,ylim=c(0,S))
}

Os argumentos desta função são:

• S: número inicial de espécies
• j: número inicial de indivíduos por espécies. Começamos com o mesmo número de indivíduos por espécie, portanto o tamanho da comunidade será $$J=Sj$$.
• D : número de mortes por ciclo, que manteremos sempre em uma.
• ciclo: número de intervalos a simular
• step: intervalo de registro dos dados, como nas funções anteriores.

Simule uma comunidades com 100 espécies e 2 indivíduos por espécie:

sim.hub1(S=100, j=2)

O que acontece com o número de espécies com o passar do tempo? Verifique se isto muda aumentando o tamanho da comunidade, que é o produto $$Sj$$. Portanto basta manter o mesmo número de espécies e aumentar o número de indivíduos por espécie:

par(mfrow=c(2,2))## para 4 gráficos na mesma janela
sim.hub1(S=100, j=2)
sim.hub1(S=100, j=4)
sim.hub1(S=100, j=8)
sim.hub1(S=100, j=12)
par(mfrow=c(1,1)) ## Volta a um grafico por janela

Sabemos que as comunidades não são sistemas fechados. Então a chegada de migrantes pode compensar a perda de espécies que vimos na simulação anterior. Vamos supor, então, que há um reservatório externo de migrantes, que chamamos metacomunidade. Uma maneira bem simples de se fazer isto é supor uma metacomunidade infinita, com todas as espécies do início da simulação, nas proporções iniciais. Precisamos definir também a taxa de migração: ela será a probabilidade de um indivíduo morto na comunidade ser substituído por um propágulo vindo de fora, da metacomunidade.

Abaixo está a função de R modificada, que inclui esta modificação. Ela tem mais um argumento, m, que é taxa de migração.

sim.hub2=function(S= 100, j=10, D=1, ciclo=2e4, step=1000, m=0.05){ 
  ## Tamanho da comunidade
  J <- S*j
  ##Matrizes para guardar os resultados
  ## matriz da especie de cada individuo por ciclo
  ind.mat=matrix(nrow=J,ncol=1+ciclo/step) 
  ##CONDICOES INICIAIS##
  ## Todas as especies comecam com o mesmo numero de individuos (j=J/S)
  ## Rotulo de especies para cada um dos inividuos
  ind.mat[,1] <- rep(1:S,each=j)
  ## Repetindo este rotulo no vetor que sofrera modificacoes
  cod.sp <- ind.mat[,1]
  ##Aqui comecam as simulacoes
  for(i in 2:(1+ciclo/step)){
    for(j in 1:step){
      ##Indice dos individuos que morrem
      morte <- sample(1:J,D)
      ## Indice dos individuos mortos que serao repostos por migrantes
      defora <- sample(c(TRUE,FALSE),size=D,replace=T,prob=c(m,1-m))
      ##Indice dos individuos que produzem filhotes para substituir os mortos
      novosd <- sample(1:J,D-sum(defora),replace=T)
      novosf <- sample(1:J,sum(defora),replace=T)
      ##Substituindo
      ## Mortos por propagulos de dentro
      if(length(novosd)>0){
        cod.sp[morte[!defora]]<-cod.sp[novosd]
      }
      ## Mortos por propagulos de fora
      if(length(novosf)>0){
        cod.sp[morte[defora]]<-ind.mat[,1][novosf]
      }
    }
    ## A cada step ciclos os resultados sao gravados
    ind.mat[,i] <- cod.sp
  }
  tempo <- seq(0,ciclo,by=step)
  colnames(ind.mat) <- tempo
  invisible(ind.mat)
  plot(tempo,apply(ind.mat,2,rich), xlab="Tempo (ciclos)", ylab="N de espécies", type="l",
       main=paste("Dinâmica Neutra com Colonização da Comunidade Original", "\n S=",S," J=",J," m=",m),ylim=c(0,S))
  }

Compare a dinâmica de número de espécies ao longo do tempo em comunidades sem migração, e com valores crescentes de taxa de migração com os comando abaixo. Em todos começamos com uma comunidade com 100 espécies, com dois indivíduos por espécies.

par(mfrow=c(2,2)) ## abre espaço para 4 graficos na mesma janela
sim.hub2(S=100, j=2,m=0)
sim.hub2(S=100, j=2,m=0.1)
sim.hub2(S=100, j=2,m=0.2)
sim.hub2(S=100, j=2,m=0.4)
par(mfrow=c(1,1))

O que acontece se aumentamos o tamanho da comunidade? Experimente começar com 10 indivíduos por espécie.

Um reservatório infinito de espécies não parece ser uma premissa muito realista. Que tal substituí-lo por um conjunto de populações com a mesma dinâmica que usamos para a comunidade? Teríamos, então, dois sistemas acoplados, cada um com sua dinâmica estocástica de nascimentos e mortes.

Mas se a metacomunidade também segue a dinâmica estocástica de soma zero, também perderá espécies com o tempo. Como resolver? Começamos por admitir que a metacomunidade é muito maior que a comunidade, pois representa o pool regional de colonizadores. Ou seja, é um sistema bem maior, pois tem mais espécies e indivíduos. Vamos supor, muito modestamente, que nela há o dobro de espécies da comunidade, cada uma com dez vezes mais indivíduos.

Apenas para lembrar o efeito do tamanho da comunidade sobre a erosão de espécies, use novamente a função de simulação sem migração para comparar sistemas que diferem nesta ordem de grandeza:

par(mfrow=c(2,1))
sim.hub1(S=100, j=2, ciclo=2e4,step=500)
sim.hub1(S=200, j=20, ciclo=2e4,step=500)
par(mfrow=c(1,1))

Já vemos que para tamanhos razoáveis (ou mesmo pequenos) de metacomunidades a erosão de espécies é bem lenta. Portanto, uma entrada de espécies também a uma taxa muito lenta já seria suficiente para compensar as extinções. Se for tão lenta quanto o tempo necessário para a evolução de uma nova espécie no sistema já temos a solução: na metacomunidade, as espécies perdidas são repostas por novas que surgem, no tempo evolutivo!

Assim, definimos uma taxa de especiação, $$nu$$, que expressa a probabilidade de um indivíduo morto na metacomunidade ser reposto por um indivíduo de uma nova espécie. Esta taxa é extremamente baixa, mas pode ser suficiente para manter, ou mesmo elevar, o número de espécies na metacomunidade.

Aqui vai a função para simular estes dois sistemas acoplados, que é o cenário imaginado por Hubbell:

sim.hub3=function(Sm=200, jm=20, S=100, j=2, m=0.01, nu=0.0001, D=1, ciclo=1e4, step=100){ 
  ## Tamanho da metacomunidade
  Jm <- Sm*jm
  ## Tamanho da comunidade
  J <- S*j
  ##Matrizes para guardar os resultados
  ## matriz da especie de cada individuo por ciclo
  ## Na metacomunidade
  meta.mat=matrix(nrow=Jm,ncol=1+ciclo/step) 
  ## Na comunidade
  ind.mat=matrix(nrow=J,ncol=1+ciclo/step)
  ##CONDICOES INICIAIS##
  ## Todas as especies comecam com o mesmo numero de individuos (j=J/S)
  ## METACOMUNIDADE
  meta.mat[,1] <- rep(1:Sm,each=jm)
  ## Repetindo este rotulo no vetor que sofrera modificacoes
  meta.sp <- meta.mat[,1]
  ##COMUNIDADE
  ## Rotulo de especies para cada um dos individuos
  ind.mat[,1] <- rep(1:S,each=j)
  ## Repetindo este rotulo no vetor que sofrera modificacoes
  cod.sp <- ind.mat[,1]
  ##Aqui comecam as simulacoes
  for(i in 2:(1+ciclo/step)){
    for(j in 1:step){
      ##Indice dos individuos que morrem
      ## Na comunidade
      morte <- sample(1:J,D)
      ## Na metacomunidade
      meta.morte <- sample(1:Jm,D)
      ## Indice dos individuos mortos da comunidade que serao repostos por migrantes 
      defora <- sample(c(TRUE,FALSE),size=D,replace=T,prob=c(m,1-m))
      ## Indice dos individuos mortos da metacomunidade que serao repostos por novas especies 
      meta.defora <- sample(c(TRUE,FALSE),size=D,replace=T,prob=c(nu,1-nu))
      ##Indice dos individuos que produzem filhotes para substituir os mortos da comunidade
      novosd <- sample(1:J,D-sum(defora),replace=T)
      novosf <- sample(1:Jm,sum(defora),replace=T)
      ##Indice dos individuos que produzem filhotes para substituir os mortos da metacomunidade
      meta.novosd <- sample(1:Jm,D-sum(meta.defora),replace=T)
      meta.novosf <- sample(1:Jm,sum(meta.defora),replace=T)
      ##Substituindo
      ## N metacomunidade ##
      ## Mortos por propagulos de dentro
      if(length(meta.novosd)>0){
        meta.sp[meta.morte[!meta.defora]]<-meta.sp[meta.novosd]
      }
      ## Mortos por novas especies
      if(length(meta.novosf)>0){
        meta.sp[meta.morte[meta.defora]]<-max(meta.sp)+1
      }
      ## Na comunidade ##
      ## Mortos por propagulos de dentro
      if(length(novosd)>0){
        cod.sp[morte[!defora]]<-cod.sp[novosd]
      }
      ## Mortos por propagulos de fora
      if(length(novosf)>0){
        cod.sp[morte[defora]]<-meta.sp[novosf]
      }
    }
    ## A cada step ciclos os resultados sao gravados
    ind.mat[,i] <- cod.sp
    meta.mat[,i] <- meta.sp
  }
  tempo <- seq(0,ciclo,by=step)
  colnames(ind.mat) <- tempo
  colnames(meta.mat) <- tempo
  resultados <- list(metacomunidade=meta.mat,comunidade=ind.mat)
  invisible(resultados)
  ## Graficos
  plot(tempo,apply(meta.mat,2,rich), xlab="Tempo (ciclos)", ylab="N de espécies", type="l",
       main=paste("Dinâmica Neutra com Colonizacao da Metacomunidade", "\n Jm=",Jm," nu=",nu," Theta=",2*Jm*nu,
         "S=",S," J=",J," m=",m), ylim=c(0,max(apply(meta.mat,2,rich))))
  lines(tempo,apply(ind.mat,2,rich),col="red")
  }

Agora temos argumentos também para os parâmetros da metacomunidade:

• Sm: número de espécies
• jm: número de indivíduos por espécie
• nu: taxa de especiação

Usando os tamanhos de comunidades e metacomunidades que já definimos, avalie o efeito de aumentar a taxa de migração, mantendo os outros parâmetros constantes:

par(mfrow=c(2,2))
sim.hub3(S=100, j=2,Sm=200,jm=20,nu=1e-9, m=0)
sim.hub3(S=100, j=2,Sm=200,jm=20,nu=1e-9, m=0.1)
sim.hub3(S=100, j=2,Sm=200,jm=20,nu=1e-9, m=0.2)
sim.hub3(S=100, j=2,Sm=200,jm=20,nu=1e-9, m=0.4)

Experimente também variar os tamanhos da comunidade e da metacomunidade, as taxas de migração e de especiação. Outra boa idéia é aumentar o tempo das simulações, para avaliar a dinâmica a longo prazo. Para isto, aumente o valor do argumento ciclo, ou a simulação pode ficar muito lenta.

• Em escala de tempo ecológico a metacomunidade desta simulação tem efeito muito diferente da metacomunidade fixa e infinita da simulação anterior?
• Qual o efeito de uma maior taxa de especiação na metacomunidade sobre a dinâmica da metacomunidade?
• O que acontece se a metacomunidade é muito pequena?

Ao invés de copiar e colar o código de cada função, você pode carregar todas elas de uma vez no R. Para isso, copie o arquivo de funções abaixo para seu diretório, e em seguida digite o comando

source("funcoes_neutr.r")

• Códigos das funções: funcoes_neutr.r
• Todos os comandos deste exercício: ex_neutra.r ³⁾