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
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.
O comportamento geral do sistema pode ser descrito pelas equações abaixo:
A primeira das equações acima pode ser interpretada da seguinte forma. Em um dado ponto, a variação na concentração
aumenta proporcionalmente ao laplaciano de
naquele ponto, i.e., quando a concentração
na vizinhança desse ponto é alta; e proporcionalmente à taxa
de reposição de
(taxa de alimentação, ou feed rate). A concentração
diminui com o termo reativo
, que representa a reação
.
De outro lado, na segunda das equações acima, a concentração
aumenta com o termo reativo
e também proporcionalmente ao laplaciano de
naquele ponto, mas diminui com a remoção de
a uma taxa
, mais rápida, portanto, do que a reposição de
.
e
são os parâmetros do modelo, juntamente com os coeficientes de difusão
e
.
O sistema pode ser imaginado como consistindo em duas substâncias
e
, envoltas por uma membrana semipermeável e imersas em um meio em que essas mesmas duas substâncias estão presentes. A membrana permite a entrada da substância
, mas não da substância
, e permite a saída da substância
, mas não da substância
.[1]
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\quad (3)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a1c3afcc002e504243e90ee4d55780147e04f780)
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\quad (4)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8b30cf254adc9bd565736d0adc73d586f6ab6a94)
Evidentemente,
é solução dessa equação, implicando em
, como já havíamos inspecionado. Alternativamente, considerando
, podemos dividir (4) 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 (3), 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.
Nesse caso, então, há três soluções estacionárias
do sistema:[2]
![{\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}})&\quad (5)\\&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/edb5faf0f9cde6fc3fac994762d20562c4df7881)
Estabilidade dos estados estacionários (sem difusão)
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}}\quad (6)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ae792c7eeb16dac2e422908aee72a26455c1c399)
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}}\quad (7)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0bc3914603b782ac3d8520aa09ba163dab45feda)
- 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\quad (8)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/477d374c9e66d4a7e233d405065b5c4f100edea6)
- 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)
- Assim, a matriz jacobiana desses pontos fica:
![{\displaystyle J_{R}(u_{i}^{*},v_{i}^{*})={\begin{bmatrix}-(v_{i}^{*})^{2}-F&-2(F+k)\\(v_{i}^{*})^{2}&(F+k)\end{bmatrix}},\quad i=1,2\quad (9)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/12bf50dd9fbbc842b1733974d9e0b74d5904c953)
- Sabemos que o produto dos autovalores dessa matriz é igual ao seu determinante. Calculando-o, obtém-se:
![{\displaystyle \Delta _{i}:=\operatorname {det} \left(J_{R}(u_{i}^{*},v_{i}^{*})\right)=(F+k)\left[(v_{i}^{*})^{2}-F\right]\quad (10)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d30ce1b2dd5a0638ce7860d81732fa18fab8be48)
- Dividindo por
:[2]
![{\displaystyle {\frac {\Delta _{i}}{F(F+k)}}={\frac {(v_{i}^{*})^{2}}{F}}-1=\left[{\frac {1\mp {\sqrt {1-4(\gamma {\sqrt {F}})^{2}}}}{2(\gamma {\sqrt {F}})}}\right]^{2}-1=\left[{\frac {1\mp {\sqrt {1-4a^{2}}}}{2a}}\right]^{2}-1}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a6dc4e89a50fb14385e69f7ae2244be3caa56cfa)
- onde se definiu
(observação: este é o
definido no Gros[2]). Nota-se que a condição de existência
para os dois pontos não-triviais é equivalente a
. Expandindo os termos, é possível mostrar que a expressão acima pode ser reescrita como:
![{\displaystyle {\frac {\Delta _{i}}{F(F+k)}}={\frac {\sqrt {1-4a^{2}}}{2a^{2}}}\left({\sqrt {1-4a^{2}}}\mp 1\right)\quad (11)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ccede7928e4bfa5bc14fcd184706e8919dc6fb97)
- Para o caso
(sinal negativo em (11)), temos a cota superior
. Portanto,
para todo
que satisfaça a condição de existência. Como o determinante é negativo, sabemos que os autovalores são reais (comentário: como as entradas da matriz são reais, se os autovalores fossem complexos, seriam também conjugados, de modo que o produto deles fosse igual ao módulo ao quadrado de qualquer um, que seria um valor positivo). Ademais, como seu produto é negativo, eles têm sinais opostos; isto é, um deles é positivo, de modo que o ponto
nunca seja estável. Depreendemos desse raciocínio que o determinante da matriz jacobiana de entradas reais ser positivo é uma condição necessária para que haja estabilidade do ponto.
- Já para
(sinal positivo em (11)), temos sempre que
. Para verificar a estabilidade, temos que agora calcular o traço da matriz jacobiana, pois o traço é a soma dos autovalores: se os autovalores são reais, eles têm o mesmo sinal por seu determinante ser positivo, de modo que o traço compartilhe o sinal com os dois autovalores; se os autovalores
são complexos, eles serão conjugados e o traço será
, de modo que a parte real dos autovalores tenha o mesmo sinal do traço. Assim, basta que o traço seja negativo para que o ponto seja estável, e que seja positivo para que seja instável.
- No caso, temos que
. Esse traço é negativo quando
e positivo quando
; ou seja,
é estável quando
e instável quando
(lembrando que
para todo
e
). Desse modo, pode-se caracterizar uma transição de estabilidade quando
.
- Utilizando simultaneamente as equações (3) e (8), obtemos:[2]
![{\displaystyle {\frac {F+k}{v_{2}^{*}}}=u_{2}^{*}=1-\gamma v_{2}^{*}=1-{\frac {F+k}{F}}v_{2}^{*}\Rightarrow (F+k){\frac {F+(v_{2}^{*})^{2}}{F}}=v_{2}^{*}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a877ef4137b90b0c22e50fc3e61db852d66e6786)
- Substituindo
, obteremos ao final:
![{\displaystyle (F+k)^{2}=F{\sqrt {k}}\quad (12)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/219c439f325d4dbb99d0ea4d3219172ce7b65d10)
Estabilidade dos estados estacionários (com difusão)
Precisamos agora analisar a estabilidade dos pontos estacionários na presença de difusão, como prescreve o sistema de equações (1), que descreve o modelo. Para isso, é necessário levar em consideração, para cada um dos estados de equilíbrio, os autovalores da matriz
, em que
é a matriz diagonal cujas entradas são
e
:[3]
![{\displaystyle D={\begin{bmatrix}D_{u}&0\\0&D_{v}\end{bmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c4d0b7564740832b583788bb1bcd2ff40f92e96c)
Se escrevermos, genericamente, que
, teremos a seguinte matriz jacobiana de reação-difusão:
![{\displaystyle \left(J_{R}-D\omega ^{2}\right){\Bigg |}_{(u_{i}^{*},v_{i}^{*})}={\begin{bmatrix}a-D_{u}\omega ^{2}&b\\c&d-D_{v}\omega ^{2}\end{bmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e679937a32845a44e2f7e93931f5e95732f7fcec)
Como já detalhado acima, para que o ponto seja estável, tal matriz tem que ter a parte real de todos os seus autovalores negativa, de modo que seu determinante seja positivo (
) e seu traço negativo (
).[4] Impondo tais condições à matriz acima, obteremos, após manipulações:[3]
![{\displaystyle {\begin{aligned}aD_{v}\omega ^{2}+dD_{u}\omega ^{2}-D_{u}D_{v}\omega ^{4}&<\operatorname {det} (J_{R})&\quad (13)\\D_{u}\omega ^{2}+D_{v}\omega ^{2}&>\operatorname {Tr} (J_{R})&\quad (14)\\\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/486e81d62f4571329e8f14ee1198a0c1592e5049)
Se o traço é negativo, vemos que a segunda equação é imediatamente satisfeita, pois o lado esquerdo é positivo em qualquer situação.[3]
- Para o ponto
, utilizamos a matriz (7), obtendo as seguintes desigualdades:
![{\displaystyle {\begin{aligned}-FD_{v}\omega ^{2}-(F+k)D_{u}\omega ^{2}-D_{u}D_{v}\omega ^{4}&<F(F+k)\\D_{u}\omega ^{2}+D_{v}\omega ^{2}&>-2F-k\\\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f7ddad0fbf6627dbbc4b014d2c3722bd00436ff2)
- Que são, evidentemente, satisfeitas, por análise simples de sinais de cada lado. Portanto, conclui-se que o ponto
é sempre estável, inclusive na presença de difusão.
- 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).[5] 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 neste modelo ocorre mesmo quando apenas o estado de equilíbrio trivial
existir.[6]
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 (ver, por exemplo, Modelo de Turing), 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
. 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
. Foram usadas condições de fronteira conforme a Figura 1.
Figura 1 - Grid para exemplificar as condições de fronteira usadas na simulação.
Cada elemento na matriz tem quatro vizinhos que são denominados por U (Up), D (Down), L (Left), R (Right). Na Figura 1, o elemento
, possui os vizinhos
,
; o elemento
possui como vizinhos
,
e
; o elemento
tem vizinhos
,
e
; e, finalmente, os vizinhos de
são
,
,
e
.
Essas condições de fronteira e a condição inicial explicada acima buscam reproduzir as mesmas condições usadas na simulação de Sayama.[7]
Resultados e discussão
As simulações abaixo reproduzem duas condições simuladas por Sayama.[7]
Simulações do Modelo de Gray-Scott para a concentração , com Falhou ao verificar gramática (erro de sintaxe): {\displaystyle (D_u, D_v) = (2\times10^{-5}, 10^{-5}). A concentração é maior nas áreas mais claras.}
|
Concentração de ![{\displaystyle u}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c3e6bb763d22c20916ed4f0bb6bd49d7470cffd8) para ![{\displaystyle (F,k)=(0.015,0.055)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b26db19fe765ecd8f27fe7d2d38f2bf743681eba) , de t=0 até t=2000.
|
Concentração de ![{\displaystyle u}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c3e6bb763d22c20916ed4f0bb6bd49d7470cffd8) para ![{\displaystyle (F,k)=(0.02,0.05)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/72e2f4ce55dc247e885e450b96a9d4be321356a7) , de t=0 até t=2000.
|
Em geral, todas as simulações apresentaram boa concordância qualitativa com os padrões simulados por
Sayama,[7] entretanto, as imagens não são isomorficamente equivalentes. Seguem abaixo alguns exemplos:
Uma explicação possível para as discrepâncias observadas é o tamanho do grid e a aplicação das condições iniciais.
Programa
Simulação do Modelo de Gray-Scott
Referências
- ↑ Reaction-Diffusion by the Gray-Scott Model: Pearson's Parametrization
- ↑ 2,0 2,1 2,2 2,3 Gros, p. 113
- ↑ 3,0 3,1 3,2 Sayama, pp. 287-289
- ↑ Sayama, p. 124
- ↑ Week 13, MCB111: Mathematics in Biology (Fall 2021)
- ↑ Gros, p. 115
- ↑ 7,0 7,1 7,2 Sayama, p. 268
Bibliografia
- C. Gros, "Complex and Adaptive Dynamical Systems". Springer-Verlag, Berlim, 2015.
- H. Sayama, "Introduction to the Modeling and Analysis of Complex Systems". Open SUNY Textbooks, Geneseo, NY, 2015.