Modelo espacialmente explícito

Anterior: Modelo de Levins aprimorado para 3 espécies | Índice: Ecologia | Próximo: -

Características

O modelo foi construído para ser coerente com o contexto ecológico da região norte da patagônia.

• ${\textstyle L\times L}$ fragmentos que pode estar destruídos, vagos ou ocupados (por uma ou mais espécies):
• ${\textstyle x_{i}\rightarrow }$ denota a fração de fragmentos ocupados por herbívoros da espécie ${\textstyle i}$ (${\textstyle 1:}$ guanaco, ${\textstyle 2:}$ ovelha);
• ${\textstyle y\rightarrow }$ fração ocupada por predadores;
• Tempo avança discretamente e os seguintes passos ocorrem de forma síncrona:
  • Colonização: um fragmento disponível pode ser colonizado pela espécie ${\textstyle \alpha \in \left\{x_{1},x_{2},y\right\}}$ por um fragmento como primeiro vizinho ocupado, com probabilidade de colonização ${\textstyle c_{\alpha }}$
  • Extinção: Um fragmento ocupado pode se tornar vago pela espécie ${\textstyle \alpha }$ por probabilidade de extinção local ${\textstyle e_{\alpha }}$
    • Predação: Um fragmento que é ocupado tanto por presa quanto um predador, tem uma probabilidade de extinção da presa, dada por uma probabilidade ${\textstyle \mu _{\alpha }}$ (${\textstyle \mu _{y}=0}$). Neste caso então é utilizado a probabilidade de extinção local ($e_{\alpha }+\mu _{\alpha }$), ou seja, consideramos como eventos mutuamente exclusivos.
  • Deslocamento competitivo: Um fragmento ocupado por ambos herbívoros, pode ser liberado do herbívoro inferior (${\textstyle x_{2}}$) com probabilidade ${\textstyle c_{x_{1}}}$.

Se temos ${\textstyle n}$ vizinhos ocupados na vizinhança de um fragmento disponível para a colonização, a probabilidade de ser colonizado por qualquer um dos vizinhos é dado por:

$$p_{\alpha }=1-\left(1-c_{\alpha }\right)^{n}$$

Uma vez que ${\textstyle b_{\alpha }=\left(1-c_{\alpha }\right)}$ é a probabilidade de que ${\textstyle \alpha }$ não colonize o fragmento. Como são eventos independentes a probabilidade de não colonizar tendo 4 vizinhos ocupados por exemplo, seria: ${\textstyle b_{\alpha }b_{\alpha }b_{\alpha }b_{\alpha }=b_{\alpha }^{4}}$, logo a probabilidade de colonizar é dada por ${\textstyle p_{\alpha }=1-b_{\alpha }^{4}}$. A dinâmica do modelo é estudada através de uma simulação computacional realizada em um sistema fechado por barreiras impenetráveis, isto é, o perímetro é composto por fragmentos destruídos, então a grade tem tamanho ${\textstyle L+2\times L+2}$. O algoritmo pode ser sintetizado nos seguintes passos:

1. Definir os parâmetros do modelo;
2. Destruir a fração ${\textstyle D}$ de fragmentos;
3. Distribuir a população inicial:
   • Cada herbívoro ocupa 50% dos fragmentos disponíveis.
   • O predador ocupa 50% dos fragmentos ocupados por herbívoros.
4. O sistema evolui de forma síncrona seguindo regras estocástica e os eventos discutidos anteriormente;
5. O sistema evolui durante um transiente (por volta de 3000 passos) até que a flutuação nas populações se estabilize;
6. Média temporal: é realizada então a média temporal (tipicamente ao longo dos últimos 2000 passos da atual simulação) da fração de fragmentos ocupados por cada espécie na atua simulação.
7. Média do ensemble: São realizadas nova simulações (tipicamente 100) com os mesmos parâmetros, e então é calculado médias do conjunto das frações ocupadas por cada espécie.

Um ensemble nada mais é que uma idealização consistindo de um largo número de cópias mentais de um sistema, onde cada uma destas cópias corresponde a um possível estado que o sistema pode estar. Basicamente uma coleção onde todas cópias possuem as mesmas propriedades, porém em estados diferentes. Dessa forma, enquanto uma média temporal (ou espacial) é a quantidade média de uma certa quantidade de um único sistema em um certo intervalo de tempo (ou espaço), a média do ensemble é quantidade média de uma certa quantidade entre muitos sistemas idênticas em um certo momento no tempo (ou espaço).

Controle de predação

Foi analisado o sistema em diferentes cenários buscando entender a dependência do estado do sistema para diferentes probabilidades de extinção de predadores e destruição do habitat. Estes parâmetros foram escolhidos pois tem origem antropogênicas. A primeira análise analisamos a a controle de predação, varia-se a probabilidade de extinção do predador (${\textstyle e_{y}}$) para obter a fração ocupada pelas três espécies após o sistema estar em equilíbrio. Os parâmetros são:

• ${\textstyle D=0.3}$
• ${\textstyle c_{x_{2}}=0.1>c_{x_{1}}=0.05}$
• ${\textstyle e_{x_{2}}=0.01<e_{x_{1}}=0.05}$
• ${\textstyle \mu _{x_{2}}=0.8\gg \mu _{x_{1}}=0.2}$
• ${\textstyle c_{y}=0.015}$
• ${\textstyle L=100}$

Além disso, para cada valor de ${\textstyle e_{y}}$ a simulação tem a duração de 5000 passos e é calculada então a média temporal de simulação para os últimos 2000 passos. É executado 100 simulações para cada valor de $e_{y}$ e então é obtido a média do ensemble.

Destruição do habitat

A segunda análise é inspirada nos três cenários discutido no contexto. cada cenário é identificado por um fração diferente de destruição do habitat:

• Baixo conflito: ${\textstyle D=0.1}$
• Médio conflito: ${\textstyle D=0.3}$
• Alto conflito: ${\textstyle D=0.5}$

Para baixo e alto conflito é executado uma simulação de 9000 passos enquanto para médio conflito é executado de 18000 passos. Porém todas compartilham do fator em comum que as ovelhas são retiradas no passo ${\textstyle 5000}$ e então é realizado a seguinte alteração de parâmetros para reiniciar a simulação:

• ${\textstyle c_{x_{1}}=0.05\rightarrow 0.1}$
• ${\textstyle e_{x_{1}}=0.05\rightarrow 0.025}$
• ${\textstyle e_{y}=0.02\rightarrow 0.015}$
• ${\textstyle \mu _{x_{1}}=0.2\rightarrow 0.3}$

Vale destacar que aqui as simulações são executadas para um valor fixo de ${\textstyle e_{y}}$, e que com o reinicio da simulação, as espécies apesar de manterem a mesma fração ocupada, acabam sendo redistribuídas aleatoriamente pelo habitat.

Resultados

Primeiro vamos comparando resultado da simulação do código disponível abaixo com os dados originais:

Comparação entre os dados originais e os gerados pelo código em Python para uma evolução temporal com $e_{y}=0.017$.

O módulo da diferença entre a fração final ocupada por cada espécie do código em Python e os dados originais foram:

• ${\textstyle \left|\Delta x_{1}\right|=3.22\times 10^{-8}}$
• ${\textstyle \left|\Delta x_{2}\right|=9.12\times 10^{-9}}$
• ${\textstyle \left|\Delta y\right|=1.15\times 10^{-9}}$

Variando então a taxa de extinção local dos pumas (${\textstyle e_{y}}$) para e calculando as médias dos ensembles temos:

A esquerda o gráfico retirado do artigo, e a direita o gerado em Python.

Porém devido à limitação de recursos computacionais, foi executado apenas 1 vez a simulação com cada valor $e_{y}$. Além disso foi reproduzido o cenário para ${\textstyle e_{y}=0.001}$.

A esquerda o gráfico retirado do artigo, e a direita o gerado em Python.

Conforme discutido, também foi simulado o cenário de conflitos:

as três imagens superiores se referem aos gráficos gerados via Python e as inferiores ás originais do artigo.

Espécie	Baixo conflito		Médio Conflito		Alto Conflito
	Artigo	Meu	Artigo	Meu	Artigo	Meu
Guanaco (${\textstyle x_{1}}$)	${\textstyle 0.35\rightarrow 0.47}$	${\textstyle 0.34\rightarrow 0.47}$	${\textstyle 0.32\rightarrow 0.45}$	${\textstyle 0.31\rightarrow 0.44}$	${\textstyle 0.03\rightarrow 0.20}$	${\textstyle 0.02\rightarrow 0.18}$
Ovelha (${\textstyle x_{2}}$)	${\textstyle 0.26\rightarrow 0.00}$	${\textstyle 0.24\rightarrow 0.00}$	${\textstyle 0.29\rightarrow 0.00}$	${\textstyle 0.29\rightarrow 0.00}$	${\textstyle 0.29\rightarrow 0.00}$	${\textstyle 0.27\rightarrow 0.00}$
Puma (${\textstyle y}$)	${\textstyle 0.25\rightarrow 0.48}$	${\textstyle 0.25\rightarrow 0.47}$	${\textstyle 0.02\rightarrow 0.17}$	${\textstyle 0.01\rightarrow 0.19}$	${\textstyle 0.00\rightarrow <0.01}$	${\textstyle 0.00\rightarrow 0.00}$

Códigos

O seguinte código em fortran foi utilizado para gerar os números pseudo-aleatórios:

! Gerador de números pseudo-aleatórios
! Jhordan Silveira de Borba
! sbjhordan@gmail.com

	program mytest1

	real RAN2  
	EXTERNAL RAN2

	IDUM= -135
	open (unit=1,file="psd.dat") 
	DO i=1,2000000000
		write(1,*),real(RAN2(IDUM))
		!print *,i
	END DO
	close(1)  

	end

	FUNCTION RAN2(IDUM)
		PARAMETER (M=714025,IA=1366,IC=150889,RM=1.4005112E-6)
		DIMENSION IR(97)
		save IR, IY
		DATA IFF /0/
		IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
			IFF=1
			IDUM=MOD(IC-IDUM,M)
			DO J=1,97
				IDUM=MOD(IA*IDUM+IC,M)
				IR(J)=IDUM
			END DO
			IDUM=MOD(IA*IDUM+IC,M)
			IY=IDUM
		ENDIF
		J=1+(97*IY)/M
		IF(J.GT.97.OR.J.LT.1)PAUSE
		IY=IR(J)
		RAN2=IY*RM
		IDUM=MOD(IA*IDUM+IC,M)
		IR(J)=IDUM
		RETURN
	END

E o seguinte código base em Python é responsável pela simulação:

# -*- coding: UTF-8 -*-

# Modelo espacialmente explícito
# Jhordan Silveira de Borba
# sbjhordan@gmail.com

#Bibliotecas
import numpy as np                                      # Biblioteca de funções matemáticas
import copy                                             # Biblioteca com funções para copiar

#CONDIÇÕES INICIAIS ------------------------------------------------------
maxt = 15000    #tiempo total de cada realizacion
Lx   = 100      #tamaño del sustrato
Ly   = 100

#Fracciones iniciales
D    = 0.3  #fraccion inicial de sitios destruidos
x10  = 0.5  #fraccion inicial de sitios ocupados por x1
x20  = 0.5  #fraccion inicial de sitios ocupados por x2
y0   = 0.5  #fraccion inicial de sitios ocupados por y

#Tasas               
cx1   = 0.05   #colonizaciones       
cx2   = 0.1        
cy    = 0.015     
ex1   = 0.05   #extinciones       
ex2   = 0.01      
ey    = 0.017
mx1y  = 0.2    #depredaciones       
mx2y  = 0.8            
depx1 = ex1+mx1y       #mortalidad aumentada=extincion+depredacion
depx2 = ex2+mx2y


g = open("psd.dat", "r") 

#Inicializaciones
s  = np.full((Lx, Ly), [1])   #matriz de habitat (1=sitio disponible, 0=destruido)
x1 = np.full((Lx, Ly), [0])   #1 = sitio ocupado, 0 sitio vazio
x2 = np.full((Lx, Ly), [0])   #1 = sitio ocupado, 0 sitio vazio
y  = np.full((Lx, Ly), [0])   #1 = sitio ocupado, 0 sitio vazio

for i in range(1,Lx-1): 
  for j in range (1,Ly-1): 
    rnd=float(g.readline())#np.random.rand()
    if (rnd < D): s[i,j]=0
    rnd=float(g.readline())#np.random.rand()
    if (rnd < x10 and s[i,j]==1): x1[i,j]=1
    rnd=float(g.readline())#np.random.rand()
    if (rnd < x20 and s[i,j]==1): x2[i,j]=1
    rnd=float(g.readline())#np.random.rand()
    if ((x1[i,j]==1 or x2[i,j]==1) and rnd < y0): y[i,j]=1
    
for i in range(Lx): #destruimos los bordes del habitat
  s[i,0]=0
  s[i,Ly-1]=0
for j in range(Ly):
  s[0,j]=0
  s[Lx-1,j]=0

#-------------------------------------------------------
#Lazo temporal
fx1=sum(sum(x1))/((Lx-2)*(Ly-2))
fx2=sum(sum(x2))/((Lx-2)*(Ly-2))
fy =sum(sum(y)) /((Lx-2)*(Ly-2))
f = open("minha.dat", "w")
f.write("	"+str(0)+"	"+str(fx1)+"	"+str(fx2)+"	"+str(fy)+"\n")

for it in range(maxt):
  if (float(it)%(float(maxt)/100)==0.):
    print("Porcentagem: "+str(it*100/maxt))  # Exibe o passo atual

  x1old=copy.copy(x1)              #poblaciones del paso anterior
  x2old=copy.copy(x2)
  yold=copy.copy(y)
  
  #Recorro la malla 
  for i in range (1,Lx-1):
    for j in range (1,Ly-1):
        
      #colonizaciones
      if (x1old[i,j]==0 and s[i,j]==1):                           #colonizacion de x1 si el sitio esta libre de x1 y no destruido:
        nvec=x1old[i-1,j]+x1old[i+1,j]+x1old[i,j-1]++x1old[i,j+1] #cantidad de vecinos ocupados
        p=1.0-(1.0-cx1)**nvec                                     #probabilidad de ser colonizado= 1 - prob de quedar libre
        rnd=float(g.readline())#np.random.rand() 
        if (rnd < p):x1[i,j]=1
      
      if (x2old[i,j]==0 and s[i,j]==1 and x1old[i,j]==0):         #colonizacion de x2 si el sitio esta libre de x1 y x2, y no destruido: 
        nvec=x2old[i-1,j]+x2old[i+1,j]+x2old[i,j-1]+x2old[i,j+1]  #cantidad de vecinos ocupados
        p=1.0-(1.0-cx2)**nvec                                     #probabilidad de ser colonizado= 1 - prob de quedar libre 
        rnd=float(g.readline())#np.random.rand()
        if (rnd < p): x2[i,j]=1                                

      ocup=0                                                       #ocupacion de herbivoros 
      if (x1old[i,j]==1 or x2old[i,j]==1): ocup=1
      if (yold[i,j]==0 and ocup==1):                               #colonizacion de y si el sitio está libre de y,y ocupado por un x 
        nvec=yold[i-1,j]+yold[i+1,j]+yold[i,j-1]+yold[i,j+1]       #cantidad de vecinos ocupados
        p=1.0-(1.0-cy)**nvec                                       #probabilidad de ser colonizado= 1 - prob de quedar libre 
        rnd=float(g.readline())#np.random.rand()
        if (rnd < p): y[i,j]=1

      #extinciones sin depredadores
      if (x1old[i,j]==1 and yold[i,j]==0):                          #extincion de x1 si el sitio esta ocupado 
        rnd=float(g.readline())#np.random.rand()
        if (rnd < ex1): x1[i,j]=0

      if (x2old[i,j]==1 and yold[i,j]==0):                          #extincion de x2 si el sitio esta ocupado 
        rnd=float(g.readline())#np.random.rand()
        if (rnd < ex2): x2[i,j]=0

      if (yold[i,j]==1):                                            #extincion de y si el sitio esta ocupado 
        rnd=float(g.readline())#np.random.rand()
        if (rnd < ey): y[i,j]=0

      #depredacion (extinción de presas donde hay depredadores)
      if (x1old[i,j]==1 and yold[i,j]==1):                          #depredacion de y sobre x1 
        rnd=float(g.readline())#np.random.rand()
        if (rnd < depx1): x1[i,j]=0

      if (x2old[i,j]==1 and yold[i,j]==1):                          #depredacion de y sobre x2 
        rnd=float(g.readline())#np.random.rand()
        if (rnd < depx2): x2[i,j]=0

      #desplazamiento por competencia (jerarquia)
      if (x1old[i,j]==1 and x2old[i,j]==1):
        rnd=float(g.readline())#np.random.rand()
        if (rnd < cx1): x2[i,j]=0
  
  fx1=sum(sum(x1))/((Lx-2)*(Ly-2))
  fx2=sum(sum(x2))/((Lx-2)*(Ly-2))
  fy =sum(sum(y)) /((Lx-2)*(Ly-2))
  f.write("	"+str(it+1)+"	"+str(fx1)+"	"+str(fx2)+"	"+str(fy)+"\n")
  
f.close
g.close() 

Além disso, foi usado o script em R para gerar os gráficos comparando dois conjunto de dados, este script é facilmente editável pra exibir o resultado de apenas uma simulação:

# Visualização de gráficos
# Jhordan Silveira de Borba
# sbjhordan@gmail.com

Dados0 <- read.table("Endereço do arquivo 0")
Dados1 <- read.table("Endereço do arquivo 1")

par(mfrow=c(1,2))
plot( c(Dados0[,1]), c(Dados0[,2]), col="black", main="Original", ylab="Fração", xlab="passo" ,type="l",ylim=c(0,0.5))
lines(c(Dados0[,1]), c(Dados0[,3]), col="blue")
lines(c(Dados0[,1]), c(Dados0[,4]), col="red")
legend("topright", c("Guanacos","Ovelhas","Pumas"), fill=c("black","blue","red"))

plot( c(Dados1[,1]), c(Dados1[,2]), col="black", main="Meu", ylab="Fração", xlab="passo" ,type="l",ylim=c(0,0.5))
lines(c(Dados1[,1]), c(Dados1[,3]), col="blue")
lines(c(Dados1[,1]), c(Dados1[,4]), col="red")
legend("topright", c("Guanacos","Ovelhas","Pumas"), fill=c("black","blue","red"))

Principal material utilizado

1. Mathematical model of livestock and wildlife: Predationand competition under environmental disturbances (Fabiana Laguna e outros, Ecological Modelling)

Referências

Anterior: Modelo de Levins aprimorado para 3 espécies | Índice: Ecologia | Próximo: -