Modelo de Gray-Scott

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 $u$ e $v$, a reação pode ser representada como

$u+2v\to 3v$

Isso significa que uma molécula da substância $u$ é transformada em uma molécula da substância $v$ por meio da ação de outras duas moléculas da substância $v$, ou seja, $v$ é 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 $u$ e $v$ mudam com o tempo e diferem em cada ponto. Por simplicidade, assume-se que a reação reversa (i.e., $3v\to u+2v$) não ocorre.

O comportamento geral do sistema pode ser descrito pelas equações abaixo:

${\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}}$

A primeira das equações acima pode ser interpretada da seguinte forma. Em um dado ponto, a variação na concentração $u$ aumenta proporcionalmente ao laplaciano de $u$ naquele ponto, i.e., quando a concentração $u$ na vizinhança desse ponto é alta; e proporcionalmente à taxa $F$ de reposição de $u$ (taxa de alimentação, ou feed rate). A concentração $u$ diminui com o termo reativo $-uv^{2}$, que representa a reação $u+2v\to 3v$.

De outro lado, na segunda das equações acima, a concentração $v$ aumenta com o termo reativo $uv^{2}$ e também proporcionalmente ao laplaciano de $v$ naquele ponto, mas diminui com a remoção de $v$ a uma taxa $F+k$, mais rápida, portanto, do que a reposição de $u$.

$F$ e $k$ são os parâmetros do modelo, juntamente com os coeficientes de difusão $D_{u}$ e $D_{v}$.

O sistema pode ser imaginado como consistindo em duas substâncias $u$ e $v$, 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 $u$, mas não da substância $v$, e permite a saída da substância $v$, mas não da substância $u$.

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 $F,k$ e dos coeficientes de difusão $D_{u},D_{v}$ 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 $(u^{*},v^{*})=(1,0)$ 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 $\partial u/\partial t=\partial v/\partial t=0$ nas equações do sistema. Fazendo isso e dispensando os termos de difusão ($D_{u}=D_{v}=0$), obtém-se o seguinte sistema de equações:

${\begin{aligned}-&uv^{2}+F(1-u)=0\\&uv^{2}-(F+k)v=0\quad (2)\\\end{aligned}}$

Somando essas duas equações, relacionamos as variáveis $u$ e $v$:

$F(1-u)-(F+k)v=0\Rightarrow u=1-\gamma v\quad (3)$

onde definiu-se o parâmetro auxiliar $\gamma ={\frac {F+k}{F}}$.

Substituindo $u$ na segunda equação do sistema (2) (e reescrevendo $F+k=\gamma F$), ficamos com:

$\left(1-\gamma v\right)v^{2}-\gamma Fv=0\Rightarrow -\gamma v^{3}+v^{2}-\gamma Fv=0\quad (4)$

Evidentemente, $v=0$ é solução dessa equação, implicando em $u=1$, como já havíamos inspecionado. Alternativamente, considerando $v\neq 0$, podemos dividir (4) por $v$, ficando com $-\gamma v^{2}+v-\gamma F=0$. Resolvendo esta equação quadrática, obtemos duas novas soluções estacionárias para $v$:

$v_{\pm }={\frac {1}{2\gamma }}(1\pm {\sqrt {1-4\gamma ^{2}F}})$

Disso, pela relação (3), temos que os valores correspondentes para $u$ são:

$u_{\mp }={\frac {1}{2}}(1\mp {\sqrt {1-4\gamma ^{2}F}})$

É 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 ( $1-4\gamma ^{2}F\geq 0$ ). Por consequência:

$4\gamma ^{2}F\leq 1\Rightarrow 4\left({\frac {F+k}{F}}\right)^{2}F\leq 1\Rightarrow F\geq 4(F+k)^{2}$, para que existam as soluções não-triviais.

Nesse caso, então, há três soluções estacionárias $(u_{i}^{*},v_{i}^{*})$ do sistema:

${\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}}$

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, $R_{i}(u,v)$. Explicitamente, analisando o sistema (1) de equações, temos que $R_{1}(u,v)=-uv^{2}+F(1-u)$ e $R_{2}(u,v)=uv^{2}-(F+k)v$. A matriz Jacobiana do sistema é então dada por:

$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)$

Analisemos a estabilidade para os três pares $(u_{i}^{*},v_{i}^{*})$ de soluções estacionárias:

• Para $(u_{0}^{*},v_{0}^{*})=(1,0)$:

$J_{R}(u_{0}^{*},v_{0}^{*})={\begin{bmatrix}-F&0\\0&-(F+k)\end{bmatrix}}\quad (7)$

Por essa ser uma matriz diagonal, os autovalores $\lambda _{i}$ são justamente as entradas das diagonais; ou seja, $\lambda _{1}=-F$ e $\lambda _{2}=-(F+k)$. Uma vez que $F$ e $k$ são parâmetros positivos, os dois autovalores são reais e negativos, e portanto o ponto $(u_{0}^{*},v_{0}^{*})$ é sempre estável.

• Para $(u_{1,2}^{*},v_{1,2}^{*})=\left({\frac {1}{2}}(1\pm {\sqrt {1-4\gamma ^{2}F}}),{\frac {1}{2\gamma }}(1\mp {\sqrt {1-4\gamma ^{2}F}})\right)$, 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 $v\neq 0$. Desse modo, se dividirmos tal equação por $v$, percebemos que ambos os pontos obedecem a:

$u_{i}^{*}v_{i}^{*}-(F+k)=0\quad (8)$

Dessa equação, podemos calcular as entradas da segunda coluna da matriz jacobiana com facilidade:

$-2u_{i}^{*}v_{i}^{*}=-2(F+k)$

$2u_{i}^{*}v_{i}^{*}-(F+k)=(F+k)$

Assim, a matriz jacobiana desses pontos fica:

$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)$

Sabemos que o produto dos autovalores dessa matriz é igual ao seu determinante. Calculando-o, obtém-se:

$\Delta _{i}:=\operatorname {det} \left(J_{R}(u_{i}^{*},v_{i}^{*})\right)=(F+k)\left[(v_{i}^{*})^{2}-F\right]\quad (10)$

Dividindo por $F(F+k)$:

${\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$

onde se definiu $a=\gamma {\sqrt {F}}$ (observação: este é o $\gamma$ definido no Gros). Nota-se que a condição de existência $a^{2}\leq 1/4$ para os dois pontos não-triviais é equivalente a $F\geq 4(F+k)^{2}$. Expandindo os termos, é possível mostrar que a expressão acima pode ser reescrita como:

${\frac {\Delta _{i}}{F(F+k)}}={\frac {\sqrt {1-4a^{2}}}{2a^{2}}}\left({\sqrt {1-4a^{2}}}\mp 1\right)\quad (11)$

• Para o caso $i=1$ (sinal negativo em (11)), temos a cota superior $1-4a^{2}<1$. Portanto, $\Delta _{1}<0$ para todo $a$ 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 $(u_{1}^{*},v_{1}^{*})$ 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 $i=2$ (sinal positivo em (11)), temos sempre que $\Delta _{2}>0$. 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 $\lambda _{i}$ são complexos, eles serão conjugados e o traço será $\operatorname {Tr} (J_{R})=2\operatorname {Re} (\lambda _{i})$, 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 $\operatorname {Tr} (J_{R}(u_{2}^{*},v_{2}^{*}))=k-(v_{2}^{*})^{2}$. Esse traço é negativo quando $(v_{2}^{*})^{2}>k$ e positivo quando $(v_{2}^{*})^{2}<k$; ou seja, $(u_{2}^{*},v_{2}^{*})$ é estável quando $v_{2}^{*}>{\sqrt {k}}$ e instável quando $v_{2}^{*}<{\sqrt {k}}$ (lembrando que $v_{2}^{*}>0$ para todo $F$ e $k$). Desse modo, pode-se caracterizar uma transição de estabilidade quando $v_{2}^{*}={\sqrt {k}}$.

Utilizando simultaneamente as equações (3) e (8), obtemos:

${\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}^{*}$

Substituindo $v_{2}^{*}={\sqrt {k}}$, obteremos ao final:

$(F+k)^{2}=F{\sqrt {k}}\quad (12)$

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 $\left(J_{R}-D\omega ^{2}\right){\Bigg |}_{(u_{i}^{*},v_{i}^{*})}$, em que $D$ é a matriz diagonal cujas entradas são $D_{u}$ e $D_{v}$:

$D={\begin{bmatrix}D_{u}&0\\0&D_{v}\end{bmatrix}}$

Se escrevermos, genericamente, que $J_{R}(u_{i}^{*},v_{i}^{*})={\begin{bmatrix}a&b\\c&d\end{bmatrix}}$, teremos a seguinte matriz jacobiana de reação-difusão:

$\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}}$

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 ($\operatorname {det} \left(J_{R}-D\omega ^{2}\right)>0$) e seu traço negativo ($\operatorname {Tr} \left(J_{R}-D\omega ^{2}\right)<0$). Impondo tais condições à matriz acima, obteremos, após manipulações:

${\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}}$

Se o traço é negativo, vemos que a segunda equação é imediatamente satisfeita, pois o lado esquerdo é positivo em qualquer situação.

• Para o ponto $(u_{0}^{*},v_{0}^{*})$, utilizamos a matriz (7), obtendo as seguintes desigualdades:

${\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}}$

Que são, evidentemente, satisfeitas, por análise simples de sinais de cada lado. Portanto, conclui-se que o ponto $(u_{0}^{*},v_{0}^{*})$ é 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). 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 $(u_{0}^{*},v_{0}^{*})=(1,0)$ existir.

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 $f(x,y,t)$:

${\frac {\partial f}{\partial t}}\approx {\frac {f(x,y,t+\Delta t)-f(x,y,t)}{\Delta t}}$

${\frac {\partial ^{2}f}{\partial x^{2}}}\approx {\frac {f(x+\Delta x,y,t)-2f(x,y,t)+f(x-\Delta x,y,t)}{\Delta x^{2}}}$

${\frac {\partial ^{2}f}{\partial y^{2}}}\approx {\frac {f(x,y+\Delta y,t)-2f(x,y,t)+f(x,y-\Delta y,t)}{\Delta y^{2}}}$

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

$\nabla ^{2}f(x,y,t)\approx {\frac {f(x+\Delta x,y,t)-2f(x,y,t)+f(x-\Delta x,y,t)}{\Delta x^{2}}}+{\frac {f(x,y+\Delta y,t)-2f(x,y,t)+f(x,y-\Delta y,t)}{\Delta y^{2}}}$

Fazendo $\Delta x=\Delta y=\Delta h$, pode-se simplificar a discretização do laplaciano para

$\nabla ^{2}f(x,y,t)={\frac {f(x+\Delta h,y,t)+f(x-\Delta h,y,t)+f(x,y+\Delta h,t)+f(x,y-\Delta h,t)-4f(x,y,t)}{\Delta h^{2}}}$

Usando a notação $f(x\pm \Delta h,y\pm \Delta h,t\pm \Delta t)\equiv f_{i\pm 1,j\pm 1}^{n\pm 1}$ é possível então escrever as equações do modelo de forma discretizada:

$u_{i,j}^{n+1}=u_{i,j}^{n}+\left[-u_{i,j}^{n}(v_{i,j}^{n})^{2}+F(1-u_{i,j}^{n})+D_{u}{\frac {u_{i+1,j}^{n}+u_{i-1,j}^{n}+u_{i,j+1}^{n}+u_{i,j-1}^{n}-4u_{i,j}^{n}}{\Delta h^{2}}}\right]\Delta t$

$v_{i,j}^{n+1}=v_{i,j}^{n}+\left[u_{i,j}^{n}(v_{i,j}^{n})^{2}-(F+k)v_{i,j}^{n}+D_{v}{\frac {v_{i+1,j}^{n}+v_{i-1,j}^{n}+v_{i,j+1}^{n}+v_{i,j-1}^{n}-4v_{i,j}^{n}}{\Delta h^{2}}}\right]\Delta t$

Utilizou-se uma rede quadrada de tamanho $69\times 69$. O estado do inicial do sistema é aquele em que todos os pontos estão no estado de equilíbrio estável trivial $(u,v)=(1,0)$, exceto o ponto central, em que é introduzida uma perturbação com $(u,v)=(0,1)$. 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 $1$, possui os vizinhos $U=D=5$, $R=L=6$; o elemento $2$ possui como vizinhos $U=1$, $R=L=9$ e $D=10$; o elemento $3$ tem vizinhos $U=D=7$, $R=1$ e $L=8$; e, finalmente, os vizinhos de $4$ são $U=3$, $D=12$, $L=11$ e $R=2$.

Essas condições de fronteira e a condição inicial explicada acima buscam reproduzir as mesmas condições usadas na simulação de Sayama.

Resultados e discussão

As simulações abaixo reproduzem duas condições simuladas por Sayama.

Simulações do Modelo de Gray-Scott para a concentração $u$, 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 $u$ para $(F,k)=(0.015,0.055)$, de t=0 até t=2000.	Concentração de $u$ para $(F,k)=(0.02,0.05)$, de t=0 até t=2000.

Em geral, todas as simulações apresentaram boa concordância qualitativa com os padrões simulados por Sayama, entretanto, as imagens não são isomorficamente equivalentes. Seguem abaixo alguns exemplos:

Figura 2

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

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.