Equação de Lotka-Volterra Competitiva Estocástica

Introdução

As Equações de Lotka-Volterra fornecem um modelo para a previsão de sistemas biológicos considerando diversas relações entre populações. Exploraremos no vigente trabalho a relação de competitividade, juntamente da adição de um termo estocástico multiplicativo. Consideraremos duas e três populações, mostrando os gráficos de evolução temporal do número de indivíduos de cada espécie e os espaços de fase. Em todas as etapas serão mostrados os resultados considerando e desconsiderando o ruído para fins comparativos. Por fim, descreveremos o método de Fokker-Planck, utilizando um exemplo simples para a demonstração.

Equações para Duas Populações

Como uma primeira análise, considerando somente duas populações distintas, vamos explorar o modelo logístico, partindo das Equações de Lotka-Volterra, de duas espécies disputando um território, que pode ser descrito pelo seguinte par de relações:

${\frac {dx_{1}}{dt}}=r_{1}x_{1}\left(1-{\frac {x_{1}+\alpha _{12}x_{2}}{K_{1}}}\right)$

${\frac {dx_{2}}{dt}}=r_{2}x_{2}\left(1-{\frac {x_{2}+\alpha _{21}x_{1}}{K_{2}}}\right)$

Nesse par, $x_{1}$ e $x_{2}$ representam as duas populações consideradas; $r_{1}$ e $r_{2}$ indicam o crescimento inerente per-capita, a capacidade de uma espécie em se reproduzir; $K_{1}$ e $K_{2}$ retratam a capacidade de carga, o número de indivíduos limite que o meio ambiente consegue suportar considerando o nicho ecológico ao que a espécie pertence; e $\alpha _{12}$ e $\alpha _{21}$, o efeito que a espécie um tem na espécie dois e vice-versa.

Para a adição de ruído ao sistema, vamos supor que o termo estocástico seja incorporado em cada equação de maneira a afetar somente a população $x_{1}$ ou a população $x_{2}$, ou seja, estamos considerando que o ruído atinja cada população de maneira proporcional a sua quantidade de indivíduos em determinado período de tempo. Como um exemplo prático, podemos pensar no termo estocástico como a presença de um agente externo que infecte uma das populações isoladamente. Claro que o efeito resultante de existir uma doença se espalhando na espécie $x_{1}$ afeta o crescimento da espécie $x_{2}$, entretanto esse efeito não é direto, sendo uma consequência da correlação entre o par de equações.

Além disso, antes de de fato colocarmos o fator estocástico nas equações, precisamos de explicar qual o tipo de ruído a ser inserido. Para tanto, será introduzido o ruído branco, que é definido a partir de suas propriedades estatísticas que configuram sua média e sua função de correlação temporal, podendo ser expressas matematicamente por

$<\xi (t)>=0$

$<\xi (t_{1})\xi (t_{2})>=\delta (t_{2}-t_{1})$

A primeira relação indica que o ruído branco possui média igual a zero, enquanto que a segunda indica o fato de os ruídos considerando diferentes tempos são independentes e, sendo assim, dá o nome branco ao ruído, pois a intensidade espectral de um processo estocástico é a Transformada de Fourier da função de correlação, sendo ela constante por se tratar de uma Delta de Dirac, logo todas as frequências estão presentes com mesma intensidade, com branco fazendo uma analogia a essa característica.

Portanto, podemos reescrever as relações supracitadas considerando o ruído como

${\frac {dx_{1}}{dt}}=r_{1}x_{1}\left(1-{\frac {x_{1}+\alpha _{12}x_{2}}{K_{1}}}\right)+x_{1}\beta {_{1}}\xi {_{1}}(t)$

${\frac {dx_{2}}{dt}}=r_{2}x_{2}\left(1-{\frac {x_{2}+\alpha _{21}x_{1}}{K_{2}}}\right)+x_{2}\beta {_{2}}\xi {_{2}}(t)$

com $\beta _{1}$ e $\beta _{2}$ sendo as amplitudes dos ruídos e $\xi _{1}(t)$ e $\xi _{2}(t)$ os ruídos brancos.

Resolução Numérica pelo Método de Itô

Para então resolvermos as equações numericamente, vamos recorrer ao cálculo diferencial estocástico. Iniciamos, para esse fim, tomando a forma de Langevin

$dX=A(X(t),t)dt+B(X(t),t)\xi (t)dt$

Deveríamos, como próxima etapa, integrar ambos os lados da equação da seguinte forma

$\int _{X}^{X+\Delta X}dX'=\int _{t}^{t+\Delta t}A(X(t'),t')dt'+\int _{t}^{t+\Delta t}B(X(t'),t')\xi (t')dt'$

Como não podemos resolvê-la pelos métodos convencionais, pois $\xi (t)$ torna a função não diferenciável apesar de ser contínua, utilizamos a relação

$B(X(t'),t')=B\left({\frac {X(t')+X(t'+\Delta t')}{2}},t'\right)$

Essa troca se justifica pelo fato de estarmos trocando o argumento de B pela média entre os valores de X no início e no fim do intervalo de integração, correspondendo a uma melhor estimação. Agora, tomando a integral, ficamos com a aproximação

$\Delta X(t)=A(X(t),t)\Delta t+B\left(X(t)+{\frac {1}{2}}\Delta X(t),t\right)\Delta W(t)$

Por seguinte, expandimos B em série de Taylor, excluímos o termo de ordem maior que dt e fazemos a substituição $\Delta t$ por $\Delta W^{2}$

$dX(t)=\underbrace {\left(A(X(t),t)+{\frac {1}{2}}{\frac {\partial B}{\partial X}}B(X(t),t)\right)} _{A^{(I)}(X(t),t)}dt+B(X(t),t)\bullet dW(t)$

Com isso, temos que a Equação diferencial de Itô é dada por

$dX(t)=A^{(I)}(X(t),t)dt+B(X(t),t)\bullet dW(t)$

Portanto, partindo da nossa relação para $x_{1}$

${\frac {dx_{1}}{dt}}=r_{1}x_{1}\left(1-{\frac {x_{1}+\alpha _{12}x_{2}}{K_{1}}}\right)+x_{1}\beta {_{1}}\xi {_{1}}(t)$

Temos, isolando $dx_{1}$

$dx_{1}=r_{1}x_{1}\left(1-{\frac {x_{1}+\alpha _{12}x_{2}}{K_{1}}}\right)dt+x_{1}\beta {_{1}}\underbrace {\xi {_{1}}(t)dt} _{dW}$

A partir disso, sabemos os coeficientes A(X(t),t) e B(X(t),t) de Langevin

$A(X_{1}(t),t)=r_{1}x_{1}\left(1-{\frac {x_{1}+\alpha _{12}x_{2}}{K_{1}}}\right)$

$B(X_{1}(t),t)=x_{1}\beta {_{1}}$

${\frac {1}{2}}{\frac {\partial B}{\partial X_{1}}}B(X_{1}(t),t)={\frac {1}{2}}\beta _{1}^{2}x_{1}$

Assim, substituindo os valores para encontrar a equação diferencial de Itô, obtemos

$dx_{1}=r_{1}x_{1}\left(1-{\frac {x_{1}+\alpha _{12}x_{2}}{K_{1}}}+{\frac {1}{2}}\beta _{1}^{2}x_{1}\right)dt+x_{1}\beta {_{1}}dW$

Analogamente, obtemos também a equação para a variável $x_{2}$

$dx_{2}=r_{2}x_{2}\left(1-{\frac {x_{2}+\alpha _{21}x_{1}}{K_{2}}}+{\frac {1}{2}}\beta _{2}^{2}x_{2}\right)dt+x_{2}\beta {_{2}}dW$

Resultados

Caso Um

Com duas populações, vamos exemplificar nesta seção o caso em que a competitividade se dá de maneira igual, ou seja, todas as constantes possuem mesmo valor independentemente da população. A diferença, então, se dará a partir do ruído, que será maior para uma delas.

Utilizando $r_{1}$ = $r_{2}$ = $\alpha _{12}$ = $\alpha _{21}$ = 1, $K_{1}$ = $K_{2}$ = 60 e população inicial de 10 indivíduos, a evolução temporal e o espaço de fases do sistema de duas populações podem ser vistos nas figuras abaixo.

Podemos verificar que, a partir da adição de um ruído maior na população 2, a população 1 cedeu quando o tempo atingiu valor próximo a 30 unidades, comportamento este que não é visto no gráfico sem ruído, no qual as duas populações convivem de maneira igual e atingindo cada uma metade de sua população limite. Realizaram-se mais de 20 simulações e todas tiveram o mesmo efeito, com a população 2 prosperando em consequência do desaparecimento da população 1.

Caso Dois

Agora, consideraremos uma pequena variação nas constantes de crescimento das populações, com um aumento de duas vezes para $r_{1}$. Todos os outros parâmetros são iguais ao Caso Um. Os resultados seguem nas figuras abaixo.

Mesmo com a população um possuindo taxa de crescimento duas vezes maior, a espécie com maior amplitude de ruído sobrevive ao final da simulação.

Equações para Três Populações

Os métodos de resolução numérica para três populações são completamente análogos aos métodos para duas populações; portanto, a matemática envolvida será omitida a partir de agora e somente serão discutidos os resultados. O conjunto de equações utilizado para três espécies distintas é

${\frac {dx_{1}}{dt}}=r_{1}x_{1}\left(1-{\frac {x_{1}+\alpha _{12}x_{2}+\alpha _{13}x_{3}}{K_{1}}}\right)+x_{1}\beta {_{1}}\xi {_{1}}(t)$

${\frac {dx_{2}}{dt}}=r_{2}x_{2}\left(1-{\frac {x_{2}+\alpha _{21}x_{1}+\alpha _{23}x_{3}}{K_{2}}}\right)+x_{2}\beta {_{2}}\xi {_{2}}(t)$

${\frac {dx_{3}}{dt}}=r_{3}x_{3}\left(1-{\frac {x_{3}+\alpha _{31}x_{1}+\alpha _{32}x_{2}}{K_{1}}}\right)+x_{3}\beta {_{3}}\xi {_{3}}(t)$

Resultados

Caso Um

Como primeira análise, vamos ao caso análogo ao anterior de duas populações. Vamos considerar espécies semelhantes, trocando entre elas somente o fator de amplitude do ruído, com o maior fator sendo da população 2. Utilizamos $K_{1}$ = $K_{2}$ = $K_{3}$ = 60, população inicial de 10 indivíduos e as outras constantes como um. As figuras a seguir mostram a evolução temporal e o espaço de fases do sistema.

Podemos verificar que, assim como na seção anterior, a população com maior fator de ruído prosperou enquanto que as outras tiveram sua extinção em um tempo de 5 unidades. Quando desconsideramos o ruído, vemos que as três espécies convivem com mesmo comportamento populacional.

Caso Dois

Para um segundo caso, vamos considerar os mesmos parâmetros de antes, modificando apenas os coeficientes entre as populações, dados na seguinte tabela.

Coeficientes $\alpha$

	$\alpha _{i1}$	$\alpha _{i2}$	$\alpha _{i3}$
$\alpha _{1j}$	1	2	1
$\alpha _{2j}$	1	1	2
$\alpha _{3j}$	2	1	1

Eles indicam que a população 1 tem dominância sobre a população 2, a população 2 tem dominância sobre a população 3 e a população 3 tem dominância sobre a população 1. As imagens podem ser vistas a seguir.

Vemos que, mesmo considerando as relações de dominância, a população com maior ruído sobressai no final, sendo a única espécie viva ao final da simulação.

Método de Fokker-Planck

A Equação de Fokker-Planck trata da equação diferencial da distribuição de probabilidade das variáveis de um sistema de Equações Diferenciais Estocásticas. Avaliaremos sua aplicação no sistema Lotka-Volterra competitivo para duas populações.

Para um processo estocástico vetorial $X=[X_{1},X_{2},...,X_{n}]$, a equação de Fokker-Planck tem a forma

${\frac {\partial P(x,t)}{\partial t}}=\sum _{i=1}^{n}{\frac {\partial [K_{i}(x,t)P(x,t)]}{\partial x_{i}}}+{\frac {1}{2}}\sum _{i,j=1}^{n}{\frac {\partial ^{2}[D_{ij}(x,t)P(x,t)]}{\partial x_{i}\partial x_{j}}}$

em que

$K_{i}(x,t)=-A_{i}(x,t)$

e

$D_{ij}(x,t)=\sum _{\alpha =1}^{m}B_{i,\alpha }(x,t)B_{j,\alpha }(x,t)$

Assim, para o nosso caso, temos

${\frac {\partial P(x_{1},x_{2},t)}{\partial t}}=-{\frac {\partial }{\partial x_{1}}}\left[r_{1}x_{1}\left(1-{\frac {x_{1}+\alpha _{12}x_{2}}{k_{1}}}\right)P(x_{1},x_{2},t)\right]-{\frac {\partial }{\partial x_{2}}}\left[r_{2}x_{2}\left(1-{\frac {x_{2}+\alpha _{21}x_{1}}{k_{2}}}\right)P(x_{1},x_{2},t)\right]+{\frac {1}{2}}\beta _{1}^{2}{\frac {\partial ^{2}x_{1}^{2}P(x_{1},x_{2},t)}{\partial x_{1}^{2}}}+{\frac {1}{2}}\beta _{2}^{2}{\frac {\partial ^{2}x_{2}^{2}P(x_{1},x_{2},t)}{\partial x_{2}^{2}}}$

Resolvendo as derivadas por regra do produto e separando nas derivadas parciais de $P(x_{1},x_{2},t)$ ficamos com

${\frac {\partial P(x_{1},x_{2},t)}{\partial t}}=\delta _{1}{\frac {\partial ^{2}P}{\partial x_{1}^{2}}}+\delta _{2}{\frac {\partial ^{2}P}{\partial x_{2}^{2}}}+\gamma _{1}{\frac {\partial P}{\partial x_{1}}}+\gamma _{2}{\frac {\partial P}{\partial x_{2}}}+\omega P$

com

$\delta _{1}={\frac {1}{2}}\beta _{1}^{2}x_{1}^{2}$

$\delta _{2}={\frac {1}{2}}\beta _{2}^{2}x_{2}^{2}$

$\gamma _{1}={\frac {r_{1}x_{1}^{2}}{k_{1}}}+{\frac {r_{1}x_{1}\alpha _{12}x_{2}}{k_{1}}}-r_{1}x_{1}+2\beta _{1}^{2}x_{1}$

$\gamma _{2}={\frac {r_{2}x_{2}^{2}}{k_{2}}}+{\frac {r_{2}x_{2}\alpha _{21}x_{1}}{k_{2}}}-r_{2}x_{2}+2\beta _{2}^{2}x_{2}$

$\omega ={\frac {2r_{1}x_{1}}{k_{1}}}+{\frac {2r_{2}x_{2}}{k_{2}}}+{\frac {r_{1}\alpha _{12}x_{2}}{k_{1}}}+{\frac {r_{2}\alpha _{21}x_{1}}{k_{2}}}-r_{1}-r_{2}+\beta _{1}^{2}+\beta _{2}^{2}$

Discretizando a equação diferencial acima, usando que

${\frac {\partial P(x,t)}{\partial t}}={\frac {P_{i,j}^{n+1}-P_{i,j}^{n}}{dt}}$

${\frac {\partial P(x,t)}{\partial x_{1}}}={\frac {1}{2}}{\frac {P_{i+1,j}^{n}-P_{i-1,j}^{n}}{dx_{1}}}$

${\frac {\partial P(x,t)}{\partial x_{2}}}={\frac {1}{2}}{\frac {P_{i,j+1}^{n}-P_{i,j-1}^{n}}{dx_{2}}}$

${\frac {\partial ^{2}P(x,t)}{\partial x_{1}^{2}}}={\frac {P_{i+1,j}^{n}+P_{i-1,j}^{n}-2P_{i,j}^{n}}{dx_{1}^{2}}}$

${\frac {\partial ^{2}P(x,t)}{\partial x_{2}^{2}}}={\frac {P_{i,j+1}^{n}+P_{i,j-1}^{n}-2P_{i,j}^{n}}{dx_{2}^{2}}}$

e que $dx_{1}=dx_{2}$, chegamos, enfim, em

$P_{i,j}^{n+1}=P_{i,j}^{n}\left(1-{\frac {2\delta _{1}dt}{dx^{2}}}-{\frac {2\delta _{2}dt}{dx^{2}}}+\omega dt\right)+P_{i+1,j}^{n}\left({\frac {\delta _{1}dt}{dx^{2}}}+{\frac {\gamma _{1}dt}{2dx}}\right)+P_{i,j+1}^{n}\left({\frac {\delta _{2}dt}{dx^{2}}}+{\frac {\gamma _{2}dt}{2dx}}\right)+P_{i-1,j}^{n}\left({\frac {\delta _{1}dt}{dx^{2}}}-{\frac {\gamma _{1}dt}{2dx}}\right)+P_{i,j-1}^{n}\left({\frac {\delta _{2}dt}{dx^{2}}}-{\frac {\gamma _{2}dt}{2dx}}\right)$

que foi resolvida numericamente.

Abaixo, temos um gráfico considerando duas populações, com constantes todas iguais e baixo ruído. O eixo das abcissas indica a população $x_{1}$, o eixo das ordenadas, a população $x_{2}$ e a cor representa a probabilidade da população possuir aquele número de indivíduos.

Código

Os principais programas utilizados no trabalho se encontram no repositório https://github.com/gugabXD/trab2-metcompC

Referências

1. Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.