Equações de Lotka-Volterra Estocásticas

Grupo: André Luis Della Valentina, Lucas dos Santos Assmann, Vinícius Bayne Müller

Introdução

O modelo de Lotka-Volterra foi desenvolvido originalmente na década de 1920, de maneira independente por Vito Volterra e Alfred Lotka, e é utilizado para descrever a dinâmica de populações com relações de predatismo. Em sua forma mais simples, as equações de Lotka-Volterra podem ser escritas como

${\begin{cases}{\dot {x}}=x(a-by)\\{\dot {y}}=y(-c+fx)\end{cases}}$

onde $x$ e $y$ denotam, respectivamente, a densidade populacional de presas e de predadores, e $a$, $b$, $c$ e $f$ são constantes positivas.

Pode-se interpretar os parâmetros da seguinte maneira:

• $a$: taxa de crescimento livre da presa;
• $b$: taxa de predação;
• $c$: taxa de mortalidade livre do predador;
• $f$: taxa de crescimento do predador devido à predação.

É interessante notar que o sistema apresenta um ponto fixo não trivial em $(x^{\ast },y^{\ast })=\left({\frac {c}{f}},{\frac {a}{b}}\right)$. Pode-se mostrar também que as demais soluções (além da trivial) são órbitas fechadas no espaço de fase.

Nesse modelo simples, não há competição entre indivíduos de uma mesma espécie e não há limite ecológico para o sustento das populações; ou seja, a população de presas cresce exponencialmente na ausência de predadores. A fim de tornar esse modelo mais realista, pretendemos estudar versões estocásticas do mesmo, que podem ser construídas de diferentes maneiras, e analisar o efeito do ruído sobre o comportamento dinâmico.

Modelo estocástico

Buscamos construir um modelo estocástico de Lotka-Volterra utilizando um processo de Wiener; contudo, a adição de ruído branco nas equações diferenciais pode ser feita de várias maneiras. Consideraremos, aqui, duas delas: ruído adicionado nos parâmetros e ruído externo ao sistema.

Ruído externo

Como primeiro modelo, simplesmente adicionaremos um ruído externo nas equações diferenciais. Esse ruído pode ser interpretado como fatores ambientais independentes das populações, causando tanto benefícios (abundância de alimento, condições climáticas favoráveis para reprodução etc.) como prejuízos (escassez, condições climáticas desfavoráveis etc.). Como as equações tratam de densidades populacionais, devemos utilizar um ruído _multiplicativo_, pois os efeitos dos fatores externos sentidos pelas populações devem ser proporcionais ao tamanho dela. Assim, escrevemos o modelo como

${\begin{cases}{\dot {X}}=X(a-bY)+\epsilon _{x}\xi _{x}X\\{\dot {Y}}=Y(-c+fX)+\epsilon _{y}\xi _{y}Y\end{cases}}$

onde os $\xi _{i}$ são os ruídos brancos de cada variável, e $\epsilon _{i}$, a intensidade dos mesmos.

Ao escrever o diferencial total, obtemos

${\begin{cases}dX=X(a-bY)dt+\epsilon _{x}X\bullet dW_{x}\\dY=Y(-c+fX)dt+\epsilon _{y}Y\bullet dW_{y}\end{cases}}$

Analisaremos o modelo proposto no sentido de Itô e no sentido de Stratonovich.

Itô

Transformando as equações acima para equações diferenciais estocásticas no sentido de Itô, obtemos

${\begin{cases}dX=X(a-bY+{\frac {1}{2}}{\epsilon _{x}}^{2})dt+\epsilon _{x}XdW_{x}\\dY=Y(-c+fX+{\frac {1}{2}}{\epsilon _{y}}^{2})dt+\epsilon _{y}YdW_{y}\\\end{cases}}$

Observa-se que, além de um incremento de Wiener proporcional a $X$ e $Y$, obtém-se parâmetros determinísticos "efetivos" $a\to a+{\frac {1}{2}}{\epsilon _{x}}^{2}$ e $c\to c-{\frac {1}{2}}{\epsilon _{y}}^{2}$.

Ruído nos parâmetros

Uma outra maneira de explorar o ruído nas equações de Lotka-Volterra é substituir o ruído de tal maneira que $a\to a+{\epsilon _{x}}X\xi _{x}(t)$ e $c\to c+{\epsilon _{y}}Y\xi _{y}(t)$, onde $\epsilon _{x}$ e $\epsilon _{y}$ são constantes e $\xi _{i}(t)$ é o ruído branco definido pelo processo de Wiener. O que modelaria um sistema que relaciona as populações entre dois níveis vizinhos da cadeia alimentar, abrangendo uma biodiversidade representada por parâmetros efetivos estocásticos nas taxas $a$ e $c$ proporcional ao número de presas e predadores, ou seja, há uma diversidade de taxas de reprodutibilidade das presas e mortalidades dos predadores, o que representaria diversas populações de especíes diferentes. Dessa maneira, obtém-se a seguinte equação diferencial:

${\begin{cases}{\dot {X}}=X(a-bY)+\epsilon _{x}\xi _{x}(t)X^{2}\\{\dot {Y}}=Y(-c+fX)-\epsilon _{y}\xi _{y}(t)Y^{2}\end{cases}}$

Seja $A_{x}(X(t),Y(t))=X(a-bY)$, $A_{y}(X(t),Y(t))=Y(-c+fX)$, $B_{x}(X(t))=\epsilon _{x}X^{2}$ e $B_{y}(Y(t))=\epsilon _{x}Y^{2}$. Integra-se essa equação em ambos os lados de tal forma que:

${\begin{cases}\int _{t}^{t+\Delta t}{dX}=\int _{t}^{t+\Delta t}{X(a-bY)dt'}+\epsilon _{x}\int _{t}^{t+\Delta t}{X^{2}\xi _{x}(t')dt'}\\\int _{t}^{t+\Delta t}{dY}=\int _{t}^{t+\Delta t}{Y(-c+fX)dt'}+\epsilon _{y}\int _{t}^{t+\Delta t}{Y^{2}\xi _{y}(t')dt'}\end{cases}}$

A parte determinista é integrada trivialmente, porém para integrar o segundo termo das equações é necessário aproximar $B_{x}(X(t))\to B_{x}(X+{\frac {1}{2}}\Delta X)$ no intervalo $\Delta t$. O que permite expandir em B em série de taylor conforme cálculo de Itô, e então obtém-se as equações diferenciais estocásticas ao tomar o limite $\Delta t\to dt$:

${\begin{cases}dX=X(a-bY+{\epsilon _{x}}^{2}X^{2})dt+\epsilon _{x}X^{2}\bullet dW_{x}\\dY=Y(-c+fX+{\epsilon _{y}}^{2}Y^{2})dt-\epsilon _{y}Y^{2}\bullet dW_{y}\\\end{cases}}$

Ao integrá-la de $t$ a $t+\Delta t$ obtém-se o método numérico utilizado para fazer as simulações computacionais. Nota-se que a integral dos segundos termos são aproximadas para $B_{i}\Delta W_{i}$ no intervalo $\Delta t$. E é utilizado o método de Heun para integrar a parte determinística da equação, obtendo por tanto o seguinte método:

${\begin{cases}X_{n+1}=X_{n}+{\frac {1}{2}}(k_{1x}+k_{2x})\Delta t+\epsilon _{x}X_{n}^{2}\Delta W_{x}\\Y_{n+1}=Y_{n}+{\frac {1}{2}}(k_{1y}+k_{2y})\Delta t-\epsilon _{y}Y_{n}^{2}\Delta W_{y}\\\end{cases}}$

Onde:

${\begin{cases}k_{1x}=X_{n}(a-bY_{n}+{\epsilon _{x}}^{2}X_{n}^{2})\\k_{2x}=(X_{n}+k_{1x}\Delta t)(a-b(Y_{n}+k_{1y}\Delta t)+{\epsilon _{x}}^{2}(X_{n}+k_{1x}\Delta t)^{2})\\k_{1y}=Y_{n}(-c+dX_{n}+{\epsilon _{y}}^{2}Y_{n}^{2})\\k_{2x}=(Y_{n}+k_{1y}\Delta t)(-c+d(X_{n}+k_{1x}\Delta t)+{\epsilon _{y}}^{2}(Y_{n}+k_{1y}\Delta t)^{2})\\\end{cases}}$

O $n$ representa o passo \Delta t na iteração do método.

Código

A estrutura do código foi a mesma para os três casos considerados, mudando apenas os detalhes da equação diferencial utilizada. A integração numérica da parte determinística foi feita pelo Método de Heun. Todas as análises foram realizadas a partir de um número N de realizações independentes e então calcula-se a média e/ou variância para cada passo de tempo destas realizações. Podendo analisar fielmente as tendências do sistema estatisticamente.

Resultados

Nesta seção, abordará-se-a diversos aspectos sobre o comportamento numérico dos sistemas simulados. Avalia-se qualitativamente estabilidade e o comportamento do sistema conforme variação dos parâmetros a, b, c, f, condições iniciais e também da intensidade do ruído.

Ruído nos parâmetros

Primeiramente realizou-se simulações variando o número de presas iniciais para um número fixo de predadores, representando-o em um mapa de fase da média de 50 realizações:

Nota-se um comportamento característico de atratores, onde o ponto fixo das equações de Lotka-Volterra determinística é o mesmo (ponto vermelho). Evidencia-se que esse ponto fixo não depende das condições iniciais.

O comportamento atrativo fica claro ao analisar a média em função do tempo:

Referências

1. BRAUER, F.; CASTILLO-CHAVEZ, C. Mathematical Models in Population Biology and Epidemiology. New York, NY: Springer New York, 2012. v. 40
2. COELHO, P. J. de O. Equações de Lotka–Volterra Estocásticas: Simulações com o Matlab. 25 de junho de 2015. Disponível em https://www.academia.edu/52185574/Sistema_de_competi%C3%A7%C3%A3o_Lotka_Volterra_sob_ru%C3%ADdo_branco. Acesso em ago. 2024.
3. SCHERER, C. Métodos computacionais da física. 2. ed ed. São Paulo: Liv. da Física, 2010.
4. KHASMINSKII, R. Z.; KLEBANER, F. C. Long term behavior of solutions of the Lotka-Volterra system under small random perturbations. The Annals of Applied Probability, v. 11, n. 3, 1 ago. 2001.