Mudanças entre as edições de "Modelo espacialmente explícito"

De Física Computacional
Ir para: navegação, pesquisa
(Criou página com '{{Ecologia| Modelo de Lotka-Volterra amortecido | - }} {{Ecologia| Modelo de Lotka-Volterra amortecido | - }}')
 
 
(9 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 1: Linha 1:
{{Ecologia| [[Modelo de Lotka-Volterra amortecido]] | - }}
+
{{Ecologia| [[Modelo de Levins aprimorado para 3 espécies]] | - }}
  
 +
=== Características ===
  
{{Ecologia| [[Modelo de Lotka-Volterra amortecido]] | - }}
+
O modelo foi construído para ser coerente com o [[contexto]] ecológico da região norte da patagônia.
 +
 
 +
*<math display="inline">L\times L</math> fragmentos que pode estar destruídos, vagos ou ocupados (por uma ou mais espécies):
 +
*<math display="inline">x_{i}\rightarrow</math> denota a fração de fragmentos ocupados por herbívoros da espécie <math display="inline">i</math> (<math display="inline">1:</math> guanaco, <math display="inline">2:</math> ovelha);
 +
*<math display="inline">y\rightarrow</math> 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 <math display="inline">\alpha\in\left\{ x_{1},x_{2},y\right\}</math> por um fragmento como primeiro vizinho ocupado, com probabilidade de colonização <math display="inline">c_{\alpha}</math>
 +
**'''Extinção''': Um fragmento ocupado pode se tornar vago pela espécie <math display="inline">\alpha</math> por probabilidade de extinção local <math display="inline">e_{\alpha}</math>
 +
***'''Predação''': Um fragmento que é ocupado tanto por presa quanto um predador, tem uma probabilidade de extinção da presa, dada por uma probabilidade <math display="inline">\mu_{\alpha}</math> (<math display="inline">\mu_{y}=0</math>). Neste caso então é utilizado a probabilidade de extinção local (<math>e_\alpha+\mu_\alpha </math>), ou seja, consideramos como eventos mutuamente exclusivos.
 +
**'''Deslocamento competitivo''': Um fragmento ocupado por ambos herbívoros, pode ser liberado do herbívoro inferior (<math display="inline">x_{2}</math>) com probabilidade <math display="inline">c_{x_{1}}</math>.
 +
 
 +
Se temos <math display="inline">n</math> 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:
 +
 
 +
<math display="block">p_{\alpha}=1-\left(1-c_{\alpha}\right)^{n}</math>Uma vez que <math display="inline">b_{\alpha}=\left(1-c_{\alpha}\right)</math> é a probabilidade de que <math display="inline">\alpha</math> não colonize o fragmento. Como são eventos independentes a probabilidade de não colonizar tendo 4 vizinhos ocupados por exemplo, seria: <math display="inline">b_{\alpha}b_{\alpha}b_{\alpha}b_{\alpha}=b_{\alpha}^{4}</math>, logo a probabilidade de colonizar é dada por <math display="inline">p_{\alpha}=1-b_{\alpha}^{4}</math>. 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 <math display="inline">L+2\times L+2</math>. O algoritmo pode ser sintetizado nos seguintes passos:
 +
#Definir os parâmetros do modelo;
 +
#Destruir a fração <math display="inline">D</math> de fragmentos;
 +
#Distribuir a população inicial:
 +
#*Cada herbívoro ocupa 50% dos fragmentos disponíveis.
 +
#*O predador ocupa 50%  dos fragmentos ocupados por herbívoros.
 +
#O sistema evolui de forma síncrona seguindo regras estocástica e os eventos discutidos anteriormente;
 +
#O sistema evolui durante um transiente (por volta de 3000 passos) até que a flutuação nas populações se estabilize;
 +
#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.
 +
#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<ref>[https://ocw.mit.edu/courses/physics/8-08-statistical-physics-ii-spring-2005/lecture-notes/microcanonical_en.pdf Statistical  Ensembles] (Instituto de Tecnologia de Massachusetts)</ref>. 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)<ref>[https://www.nii.ac.jp/qis/first-quantum/e/forStudents/lecture/pdf/noise/chapter1.pdf Mathematical Methods] (Yoshihisa Yamamoto, Instituto Nacional de Informática)</ref>.
 +
 
 +
====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 (<math display="inline">e_{y}</math>) para obter a fração ocupada pelas três espécies após o sistema estar em equilíbrio. Os parâmetros são:
 +
 
 +
*<math display="inline">D=0.3</math>
 +
*<math display="inline">c_{x_{2}}=0.1>c_{x_{1}}=0.05</math>
 +
*<math display="inline">e_{x_{2}}=0.01<e_{x_{1}}=0.05</math>
 +
*<math display="inline">\mu_{x_{2}}=0.8\gg\mu_{x_{1}}=0.2</math>
 +
*<math display="inline">c_{y}=0.015</math>
 +
*<math display="inline">L=100</math>
 +
 
 +
Além disso, para cada valor de <math display="inline">e_{y}</math> 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 <math>e_y</math> 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: <math display="inline">D=0.1</math>
 +
*Médio conflito: <math display="inline">D=0.3</math>
 +
*Alto conflito: <math display="inline">D=0.5</math>
 +
 
 +
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 <math display="inline">5000</math> e então é realizado a seguinte alteração de parâmetros para reiniciar a simulação:
 +
 
 +
*<math display="inline">c_{x_{1}}=0.05\rightarrow0.1</math>
 +
*<math display="inline">e_{x_{1}}=0.05\rightarrow0.025</math>
 +
*<math display="inline">e_{y}=0.02\rightarrow0.015</math>
 +
*<math display="inline">\mu_{x_{1}}=0.2\rightarrow0.3</math>
 +
 
 +
Vale destacar que aqui as simulações são executadas para um valor fixo de <math display="inline">e_{y}</math>, 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:
 +
 
 +
<div class="center">[[Ficheiro:comparacao.png|centro|miniaturadaimagem|Comparação entre os dados originais e os gerados pelo código em Python para uma evolução temporal com <math>e_y=0.017</math>.|alt=|646x314px]]
 +
</div>
 +
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:
 +
 
 +
*<math display="inline">\left|\Delta x_{1}\right|=3.22\times10^{-8}</math>
 +
*<math display="inline">\left|\Delta x_{2}\right|=9.12\times10^{-9}</math>
 +
*<math display="inline">\left|\Delta y\right|=1.15\times10^{-9}</math>
 +
 
 +
Variando então a taxa de extinção local dos pumas (<math display="inline">e_{y}</math>) para e calculando as médias dos ensembles temos:
 +
[[Ficheiro:media_esemble.png|centro|miniaturadaimagem|729x318px|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 <math>e_y</math>. Além disso foi reproduzido o cenário para <math display="inline">e_{y}=0.001</math>.
 +
[[Ficheiro:0001.png|centro|miniaturadaimagem|614x300px|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:<div class="center">[[Ficheiro:conflitos.png|centro|miniaturadaimagem|397x300px|as três imagens superiores se referem aos gráficos gerados via Python e as inferiores ás originais do artigo.]]
 +
</div>
 +
<div align="center">
 +
 
 +
{| class="wikitable"
 +
| align="center" |Espécie
 +
| align="center" |Baixo conflito
 +
| align="center" |
 +
| align="center" |Médio Conflito
 +
| align="center" |
 +
| align="center" | Alto Conflito
 +
| align="center" |
 +
|-
 +
| align="center" |
 +
| align="center" |Artigo
 +
| align="center" |Meu
 +
| align="center" |Artigo
 +
| align="center" | Meu
 +
| align="center" | Artigo
 +
| align="center" |Meu
 +
|-
 +
| align="center" |Guanaco (<math display="inline">x_{1}</math>)
 +
| align="center" |<math display="inline">0.35\rightarrow0.47</math>
 +
| align="center" |<math display="inline">0.34\rightarrow0.47</math>
 +
| align="center" |<math display="inline">0.32\rightarrow0.45</math>
 +
| align="center" |<math display="inline">0.31\rightarrow0.44</math>
 +
| align="center" |<math display="inline">0.03\rightarrow0.20</math>
 +
| align="center" |<math display="inline">0.02\rightarrow0.18</math>
 +
|-
 +
| align="center" | Ovelha (<math display="inline">x_{2}</math>)
 +
| align="center" |<math display="inline">0.26\rightarrow0.00</math>
 +
| align="center" |<math display="inline">0.24\rightarrow0.00</math>
 +
| align="center" |<math display="inline">0.29\rightarrow0.00</math>
 +
| align="center" |<math display="inline">0.29\rightarrow0.00</math>
 +
| align="center" |<math display="inline">0.29\rightarrow0.00</math>
 +
| align="center" |<math display="inline">0.27\rightarrow0.00</math>
 +
|-
 +
| align="center" |Puma (<math display="inline">y</math>)
 +
| align="center" |<math display="inline">0.25\rightarrow0.48</math>
 +
| align="center" |<math display="inline">0.25\rightarrow0.47</math>
 +
| align="center" |<math display="inline">0.02\rightarrow0.17</math>
 +
| align="center" |<math display="inline">0.01\rightarrow0.19</math>
 +
| align="center" |<math display="inline">0.00\rightarrow<0.01</math>
 +
| align="center" |<math display="inline">0.00\rightarrow0.00</math>
 +
|}
 +
 
 +
</div>
 +
=== Códigos===
 +
O seguinte código em fortran foi utilizado para gerar os números pseudo-aleatórios:<pre>
 +
! 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
 +
</pre>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:<pre>
 +
# 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"))
 +
</pre>
 +
 
 +
===Principal material utilizado===
 +
 
 +
#[https://arxiv.org/pdf/1409.0024.pdf Mathematical model of livestock and wildlife: Predationand competition under environmental disturbances] (Fabiana Laguna e outros, Ecological Modelling)
 +
 
 +
===Referências===
 +
<references />
 +
{{Ecologia| [[Modelo de Levins aprimorado para 3 espécies]] | - }}

Edição atual tal como às 19h45min de 13 de abril de 2021

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