Modelo espacialmente explícito

De Física Computacional
Ir para navegação Ir para pesquisar

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.

  • fragmentos que pode estar destruídos, vagos ou ocupados (por uma ou mais espécies):
  • denota a fração de fragmentos ocupados por herbívoros da espécie ( guanaco, ovelha);
  • 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 por um fragmento como primeiro vizinho ocupado, com probabilidade de colonização
    • Extinção: Um fragmento ocupado pode se tornar vago pela espécie por probabilidade de extinção local
      • Predação: Um fragmento que é ocupado tanto por presa quanto um predador, tem uma probabilidade de extinção da presa, dada por uma probabilidade (). Neste caso então é utilizado a probabilidade de extinção local (), ou seja, consideramos como eventos mutuamente exclusivos.
    • Deslocamento competitivo: Um fragmento ocupado por ambos herbívoros, pode ser liberado do herbívoro inferior () com probabilidade .

Se temos 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:

Uma vez que é a probabilidade de que não colonize o fragmento. Como são eventos independentes a probabilidade de não colonizar tendo 4 vizinhos ocupados por exemplo, seria: , logo a probabilidade de colonizar é dada por . 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 . O algoritmo pode ser sintetizado nos seguintes passos:

  1. Definir os parâmetros do modelo;
  2. Destruir a fração 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[1]. 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)[2].

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 () para obter a fração ocupada pelas três espécies após o sistema estar em equilíbrio. Os parâmetros são:

Além disso, para cada valor de 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 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:
  • Médio conflito:
  • Alto conflito:

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 e então é realizado a seguinte alteração de parâmetros para reiniciar a simulação:

Vale destacar que aqui as simulações são executadas para um valor fixo de , 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 .

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:

Variando então a taxa de extinção local dos pumas () 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 . Além disso foi reproduzido o cenário para .

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 ()
Ovelha ()
Puma ()

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

  1. Statistical Ensembles (Instituto de Tecnologia de Massachusetts)
  2. Mathematical Methods (Yoshihisa Yamamoto, Instituto Nacional de Informática)

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