Introdução
Descrição do Modelo
O modelo de Gray-Scott descreve uma reação autocatalítica. Sejam duas substâncias químicas cujas concentrações em um dado ponto do espaço são dadas pelas variáveis
e
, a reação pode ser representada como
![{\displaystyle u+2v\to 3v}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ec1c9bbcb954b57f08fa68235b4ebda25c1a26bd)
Isso significa que uma molécula da substância
é transformada em uma molécula da substância
por meio da ação de outras duas moléculas da substância
, ou seja,
é um catalisador de sua própria produção (daí o termo autocatálise). Além dessa reação, ambas substâncias se difundem pelo meio (por isso esse modelo pertence à classe mais geral de modelos reativos-difusivos) e, portanto, as concentrações
e
mudam com o tempo e diferem em cada ponto. Por simplicidade, assume-se que a reação reversa (i.e.,
) não ocorre. Há reposição de
a uma taxa
(taxa de alimentação, feed rate) e remoção de
a uma taxa ligeiramente mais rápida do que a reposição de
.
O comportamento geral do sistema pode ser descrito pelas equações abaixo:
![{\displaystyle {\begin{aligned}{\frac {\partial {u}}{\partial {t}}}&=-uv^{2}+F(1-u)+D_{u}\nabla ^{2}u\\{\frac {\partial {v}}{\partial {t}}}&=uv^{2}-(F+k)v+D_{v}\nabla ^{2}v\quad (1)\\\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f1d6d5474c9e3e54de79afe9bf25dad47886cefd)
Análise de estabilidade
Nota: A análise em toda esta seção pressupõe sempre que os parâmetros e coeficientes de difusão são positivos.
Soluções estacionárias sem difusão
O modelo de Gray-Scott depende dos parâmetros
e dos coeficientes de difusão
das espécies químicas. Ignorando em um primeiro momento os termos de difusão, percebe-se que, por inspeção, o sistema possui uma solução estacionária em
para quaisquer valores dos parâmetros. Esse ponto, no entanto, não é a única solução estacionária do sistema; para encontrar as outras, é necessário impor
nas equações do sistema. Fazendo isso e dispensando os termos de difusão (
), obtém-se o seguinte sistema de equações:
![{\displaystyle {\begin{aligned}-&uv^{2}+F(1-u)=0\\&uv^{2}-(F+k)v=0\quad (2)\\\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7b11de91d0874cfcde3a857b086e88f857e34838)
Somando essas duas equações, relacionamos as variáveis
e
:
![{\displaystyle F(1-u)-(F+k)v=0\Rightarrow u=1-\gamma v}](https://wikimedia.org/api/rest_v1/media/math/render/svg/79e84a3598f6173e5258260920d549acbca0610b)
onde definiu-se o parâmetro auxiliar
.
Substituindo
na segunda equação do sistema (2) (e reescrevendo
), ficamos com:
![{\displaystyle \left(1-\gamma v\right)v^{2}-\gamma Fv=0\Rightarrow -\gamma v^{3}+v^{2}-\gamma Fv=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c0b43014a275d7a46ccecde177526bb9d09f4ceb)
Evidentemente,
é solução dessa equação, implicando em
, como já havíamos inspecionado. Alternativamente, considerando
, podemos dividir a expressão acima por
, ficando com
. Resolvendo esta equação quadrática, obtemos duas novas soluções estacionárias para
:
![{\displaystyle v_{\pm }={\frac {1}{2\gamma }}(1\pm {\sqrt {1-4\gamma ^{2}F}})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/06ac57f2e011f42220fa2aff4dcb770028e3f237)
Disso, pela relação
, temos que os valores correspondentes para
são:
![{\displaystyle u_{\mp }={\frac {1}{2}}(1\mp {\sqrt {1-4\gamma ^{2}F}})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f071866f2e40aed74a63563ce20347da9ec74226)
É necessário apontar que, para que as duas últimas soluções (não-triviais) existam — isto é, sejam números reais — o fator dentro da raiz quadrada tem de ser positivo (
). Por consequência:
, para que existam as soluções não-triviais.
Portanto, há três soluções estacionárias
do sistema:[1]
![{\displaystyle {\begin{aligned}&u_{0}^{*}=1&v_{0}^{*}=0\\&u_{1}^{*}={\frac {1}{2}}(1+{\sqrt {1-4\gamma ^{2}F}})&v_{1}^{*}={\frac {1}{2\gamma }}(1-{\sqrt {1-4\gamma ^{2}F}})\\&u_{2}^{*}={\frac {1}{2}}(1-{\sqrt {1-4\gamma ^{2}F}})&v_{2}^{*}={\frac {1}{2\gamma }}(1+{\sqrt {1-4\gamma ^{2}F}})\\\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c7d9419a34823054dba14aae9b66920f56f08547)
Estabilidade dos estados estacionários
Para avaliar a estabilidade das soluções acima, faz-se necessário obter a matriz Jacobiana dos termos de reação,
. Explicitamente, analisando o sistema (1) de equações, temos que
e
. A matriz Jacobiana do sistema é então dada por:
![{\displaystyle J_{R}(u,v)={\begin{bmatrix}\partial _{u}R_{1}&\partial _{v}R_{1}\\\partial _{u}R_{2}&\partial _{v}R_{2}\end{bmatrix}}={\begin{bmatrix}-v^{2}-F&-2uv\\v^{2}&2uv-(F+k)\end{bmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b69a5983ff66647bf1d19c85262f3d112dd02f68)
Analisemos a estabilidade para os três pares
de soluções estacionárias:
- Para
:
![{\displaystyle J_{R}(u_{0}^{*},v_{0}^{*})={\begin{bmatrix}-F&0\\0&-(F+k)\end{bmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/067ba08696c52ddc306666b4b197a7373df14775)
- Por essa ser uma matriz diagonal, os autovalores
são justamente as entradas das diagonais; ou seja,
e
. Uma vez que
e
são parâmetros positivos, os dois autovalores são reais e negativos, e portanto o ponto
é sempre estável.
- Para
, podemos utilizar uma estratégia que simplifica as contas. Em particular, nota-se que os dois pontos obedecem à segunda equação do sistema (2) com
. Desse modo, se dividirmos tal equação por
, percebemos que ambos os pontos obedecem a:
![{\displaystyle u_{i}^{*}v_{i}^{*}-(F+k)=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6a585eec4fbe7c3d30af86014ff779898af68c9f)
- Dessa equação, podemos calcular as entradas da segunda coluna da matriz jacobiana com facilidade:
![{\displaystyle -2u_{i}^{*}v_{i}^{*}=-2(F+k)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/550e18bb317488d949cfa911a2af0cb249adb60f)
![{\displaystyle 2u_{i}^{*}v_{i}^{*}-(F+k)=(F+k)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d54f976c7beccf7fbc4d4e7dcbef16fe2b3492d4)
- Já para a primeira coluna, podemos utilizar a equação quadrática a qual
e
obedecem,
. Dela, por manipulações elementares, obtemos uma forma mais simples de calcular as entradas restantes da matriz:
![{\displaystyle -(v_{i}^{*})^{2}-F=-{\frac {v_{i}^{*}}{\gamma }}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b032755b388c647993383a3b6d3b591e06d98f90)
![{\displaystyle -(v_{i}^{*})^{2}=-{\frac {v_{i}^{*}}{\gamma }}+F}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ad9a2e0bbd75c594dcd01d5cc558f6fb6e047d2b)
- Assim, a matriz jacobiana desses pontos fica:
![{\displaystyle J_{R}(u_{i}^{*},v_{i}^{*})={\begin{bmatrix}-{\frac {v_{i}^{*}}{\gamma }}&-2(F+k)\\-{\frac {v_{i}^{*}}{\gamma }}+F&(F+k)\end{bmatrix}},\quad i=1,2.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/28d1df098fc3a9499140aac1f6a3d8b9a94a5683)
Esse estado de equilíbrio é estável porque a matriz jacobiana
possui traço negativo e determinante positivo[2].
Se agora incluímos os termos de difusão
e
, deve-se levar em consideração a matriz
. Aqui,
é a matriz jacobiana dos termos de reação,
é a matriz diagonal dos termos de difusão e
é o parâmetro que determina a frequência espacial das perturbações. A demonstração da validade desse método pode ser encontrada na referência[2]. Aplicando ao modelo de Gray-Scott em
:
Para que o estado de equilíbrio
seja estável é necessário que o determinante da matriz acima seja positivo e o traço seja negativo. Obtém-se então
Ambas desigualdades são imediatamente satisfeitas para quaisquer valores de
, e
. Portanto, o estado de equilíbrio
permanece estável no modelo de Gray-Scott mesmo após a inclusão dos coeficientes de difusão, sejam quais forem os valores desses coeficientes (lembrando que estamos nos restringindo a valores positivos dos parâmetros e coeficientes).
Esse é um resultado à primeira vista surpreendente. Em geral, o surgimento de padrões complexos e não homogêneos em sistemas reativos-difusivos está relacionado à desestabilização de um ou mais estados de equilíbrio homogêneo causada pela introdução dos coeficientes de difusão (conhecida como instabilidade de Turing)[3].
Entretanto, no caso do modelo de Gray-Scott, o surgimento de padrões complexos e não homogêneos não decorre da instabilidade de Turing, uma vez que o surgimento de padrões não triviais nesse modelo ocorre mesmo quando apenas o estado de equilíbrio trivial
está presente [1].
Implementação
Será usado o método FTCS (Foward Time Central Space) para integrar as equações do modelo. Como existem explicações do método em toda literatura e em outras entradas da Wiki Título do link, a explicação aqui será sucinta.
O método consiste em discretizar a derivada parcial em relação ao tempo para frente e discretizar as derivadas parciais de segunda ordem em relação ao espaço centralmente. Para uma função
:
A partir das duas últimas equações acima é fácil mostrar que o laplaciano em duas dimensões, como será usado no presente trabalho, pode ser escrito como
Fazendo
, pode-se simplificar a discretização do laplaciano para
Usando a notação
é possível então escrever as equações do modelo de forma discretizada:
Utilizou-se uma rede quadrada de tamanho
com condições de contorno periódicas. O estado do inicial do sistema é aquele em que todos os pontos estão no estado de equilíbrio estável trivial
, exceto o ponto central, em que é introduzida uma perturbação com
, como simulado por Sayama[2].
A formação de padrões no modelo depende fortemente não apenas dos parâmetros e coeficientes de difusão, mas também da resolução,
.
Referências
- ↑ 1,0 1,1 C. Gros, "Complex and Adaptive Dynamical Systems". Springer-Verlag, Berlim, 2015.
- ↑ 2,0 2,1 2,2 H. Sayama, "Introduction to the Modeling and Analysis of Complex Systems". Open SUNY Textbooks, Geneseo, NY, 2015.
- ↑ Week 13, MCB111: Mathematics in Biology (Fall 2021)