Modelo Brusselator de Reação-Difusão

Grupo: Carolina Lenzi, Eric Naiber e Vitória Xavier

O objetivo deste trabalho é implementar o modelo de reação-difusão Brusselator em duas dimensões, frequentemente utilizado para estudar sistemas complexos químicos e biológicos. O modelo é um sistema não linear de equações diferenciais parciais e foi proposto em 1968 por Ilya Prigogine e seus colaboradores da Universidade Livre de Bruxelas. Desde então tem sido aplicado para analisar reações oscilatórias e autocatalíticas. O método computacional utilizado para implementar o modelo foi o método FTCS (Forward Time Centered Space).

Modelo de Brusselator

O estudo de sistemas químicos e biológicos frequentemente requer o uso de modelos que caracterizam reações de reação-difusão. Um dos modelos mais utilizados é o modelo de Brusselator, que é utilizado para descrever o mecanismo químico de reação-difusão com oscilações não lineares. . Turing observou que quando determinadas reações são associadas a difusão, é possível obter um padrão espacial estável, e isso leva a teoria de morfogênese. Além de processos de reação-difusão, o modelo Brusselator é observado em reações enzimáticas e na física de plasma e de lasers.

O mecanismo de Brusselator proposto por Prigogine é dado por:

$A\to U$ (1.a)

$B+U\to V+D$ (1. b)

$2U+V\to 3U$ (1.c)

$U\to E$ (1.d)

Onde U e V são as espécies químicas de interesse. Assumimos A e B em excesso para que o sistema não atinja o estado de equilíbrio. Esse sistema químico foi importante para o avanço na área de sistemas complexos porque possibilita o uso de modelos matemáticos de duas dimensões, já que U e V são variáveis dependentes, e admite oscilações de ciclo limite.

As equações diferenciais parciais associadas com o sistema Brusselator são dadas por:

$\frac{\partial u}{\partial t}=a+u^{2}v-(b+1)u+D_u(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2})$

$\frac{\partial v}{\partial t}=bu-u^{2}v+D_v(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2})$

onde $u(x,y,t)$ e $v(x,y,t)$ são as concentrações a serem investigadas em função de tempo e espaço, $a$ e $b$ são constantes relativas às concentrações dos reagentes A e B, e $D_u$ e $D_v$ são as constantes de difusão de U e V.

A solução analítica do sistema reação-difusão Brusselator ainda não é conhecida e por isso há o interesse de explorá-la numericamente.

Análise da estabilidade do sistema

Análise de ponto crítico

Considerando o sistema livre de difusão, quando $D_u=D_v=0$:

$\frac{du}{dt}=f(u,v)=a+u^{2}v-(b+1)u$

$\frac{dv}{dt}=g(u,v)=bu-u^{2}v$

Encontramos os pontos críticos do sistema igualando as derivadas à zero. Obtemos que o único ponto crítico é $(u^{*},v^{*})=(a,\frac{b}{a})$. A matriz jacobiana $J^{*}$ no ponto crítico $(u^{*},v^{*})$ é dada por

$J^{*}=\left[\begin{array}{cc}\frac{\partial f}{\partial u}{|}_{u^{*},v^{*}}& \frac{\partial f}{\partial v}{|}_{u^{*},v^{*}}\\ \frac{\partial g}{\partial u}{|}_{u^{*},v^{*}}& \frac{\partial g}{\partial v}{|}_{u^{*},v^{*}}\end{array}\right]=\left[\begin{array}{cc}b-1& a^2\\ -b& -a^2\end{array}\right]$

Os autovalores de $J^{*}$ são os valores $\lambda$ que satisfazem a equação caracterísitca

${\lambda }^2+(1-b+a^2)\lambda +a^2=0$

Os autovalores claramente mostram dependência em $1-b+a^2$ e no determinante $\mathrm{\Delta }=(1-b+a^2{)}^2-4a^2$. Esses autovalores governam a estabilidade do ponto crítico e determinam a existência de um ciclo limite. Conforme discutido em Twizell, quando $Re(\lambda )<0$ o sistema é estável e não existe ciclo limite. Já para $Re(\lambda )>0$ o sistema é instável e existe ciclo limite.

Análise de ponto fixo

Ainda em Twizell, foi feita a análise de ponto fixo. O único ponto fixo encontrado é $(u^{*},v^{*})=(a,\frac{b}{a})$. Manipulando a matriz Jacobiana nesse ponto e analisando os autovalores, obtém-se que a solução é atraída para o ponto $(a,\frac{b}{a})$ se $1-b+a^2>0$. Caso contrário, a solução não converge.

Método FTCS

O FTCS (Forward Time Centered Space) é um método de diferença finita que utiliza a derivada à direita ("para frente") no tempo e a derivada segunda centralizada no espaço para discretizar as variáveis. As derivadas no tempo e no espaço bidimensional ficam:

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

$\frac{{\partial }^{2}f(x,y,t)}{\partial x^2}\to \frac{f(x-\mathrm{\Delta }x,y,t)-2f(x,y,t)+f(x+\mathrm{\Delta }x,y,t)}{(\mathrm{\Delta }x{)}^2}+{𝒪}(\mathrm{\Delta }{x}^3)$

$\frac{{\partial }^{2}f(x,y,t)}{\partial y^2}\to \frac{f(x,y-\mathrm{\Delta }y,t)-2f(x,y,t)+f(x,y+\mathrm{\Delta }y,t)}{(\mathrm{\Delta }y{)}^2}+{𝒪}(\mathrm{\Delta }{y}^3)$

Substituindo nas equações do Brusselator

$\frac{u(x,y,t+\mathrm{\Delta }t)-u(x,y,t)}{\mathrm{\Delta }t}=f(u,v)+D_u(\frac{u(x-\mathrm{\Delta }x,y,t)-2u(x,y,t)+u(x+\mathrm{\Delta }x,y,t)}{(\mathrm{\Delta }x{)}^2}+\frac{u(x,y-\mathrm{\Delta }y,t)-2u(x,y,t)+u(x,y+\mathrm{\Delta }y,t)}{(\mathrm{\Delta }y{)}^2})$

$\frac{v(x,y,t+\mathrm{\Delta }t)-v(x,y,t)}{\mathrm{\Delta }t}=g(u,v)+D_v(\frac{v(x-\mathrm{\Delta }x,y,t)-2v(x,y,t)+v(x+\mathrm{\Delta }x,y,t)}{(\mathrm{\Delta }x{)}^2}+\frac{v(x,y-\mathrm{\Delta }y,t)-2v(x,y,t)+v(x,y+\mathrm{\Delta }y,t)}{(\mathrm{\Delta }y{)}^2})$

onde $f(u,v)$ e $g(u,v)$ são as funções que representam a reação sem difusão.

Utilizamos discretização do tipo

$t=0,1\mathrm{\Delta }t,2\mathrm{\Delta }t,3\mathrm{\Delta }t,...,t_{max}$

$x=0,1\mathrm{\Delta }x,2\mathrm{\Delta }x,3\mathrm{\Delta }x,...,N_x$

$y=0,1\mathrm{\Delta }y,2\mathrm{\Delta }y,3\mathrm{\Delta }y,...,N_y$

Utilizando a notação $f(i\mathrm{\Delta }x,j\mathrm{\Delta }y,n\mathrm{\Delta }t)=f_{i,j}^n$, assumindo $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }s$ e rearranjando os termos, reescrevemos as equações como

$u_{i,j}^{n+1}=u_{i,j}^n+f(u_{i,j}^n,v_{i,j}^n)\mathrm{\Delta }t+K_u(u_{i-1,j}^n+u_{i+1,j}^n+u_{i,j-1}^n+u_{i,j+1}^n-4u_{i,j}^n)$

$v_{i,j}^{n+1}=v_{i,j}^n+g(u_{i,j}^n,v_{i,j}^n)\mathrm{\Delta }t+K_v(v_{i-1,j}^n+v_{i+1,j}^n+v_{i,j-1}^n+v_{i,j+1}^n-4v_{i,j}^n)$

onde $K_u=\frac{D_u\mathrm{\Delta }t}{(\mathrm{\Delta }s{)}^2}$ e $K_v=\frac{D_v\mathrm{\Delta }t}{(\mathrm{\Delta }s{)}^2}$.

Análise de estabilidade

A análise de estabilidade do método FTCS pode ser feita com a análise local de von Neumann. Para isso precisamos linearizar as equações do modelo Brusselator e escrever as soluções em modos de Fourier, da seguinte forma:

$\delta u_{i,j}^{n+1}=\delta u_{i,j}^n+K_u(\delta u_{i+1,j}^n+\delta u_{i-1,j}^n+\delta u_{i,j+1}^n+\delta u_{i,j-1}^n-4\delta u_{i,j}^n)+\mathrm{\Delta }t{f}_u\delta u_{i,j}^n+\mathrm{\Delta }t{f}_v\delta v_{i,j}^n$

$\delta v_{i,j}^{n+1}=\delta v_{i,j}^n+K_v(\delta v_{i+1,j}^n+\delta v_{i-1,j}^n+\delta v_{i,j+1}^n+\delta v_{i,j-1}^n-4\delta v_{i,j}^n)+\mathrm{\Delta }t{g}_u\delta u_{i,j}^n+\mathrm{\Delta }t{g}_v\delta v_{i,j}^n$

$(\genfrac{}{}{0ex}{}{\delta u^n}{\delta v^n})={\xi }^n{e}^{i(k_{x}i\mathrm{\Delta }x+k_{y}j\mathrm{\Delta }y)}(\genfrac{}{}{0ex}{}{\delta u_0}{\delta v_0})$

Para que o método seja estável, é preciso que $|\xi |<1$. Após inserir as soluções nas equações linearizadas e realizar manipulações algébricas, como feito por Scholz, obtemos

${\xi }_{\pm }=\frac{1}{2}[\beta \pm \sqrt{{\beta }^2-4(\mathrm{\Delta }t{)}^2{a}^{2}b}]$

onde

$\beta =2-4\gamma (K_u+K_v)+\mathrm{\Delta }t(b-1-a^2)$

e

$\gamma =si{n}^2\left(\frac{k_x\mathrm{\Delta }s}{2}\right)+si{n}^2\left(\frac{k_y\mathrm{\Delta }s}{2}\right)$

Os valores $\mathrm{\Delta }t=0.01$ e $\mathrm{\Delta }s=1$ satisfazem a condição de estabilidade, como mostrado em Scholz.

Implementação

A seguir, o trecho do código relativo à implementação do método na linguagem Python. O código completo encontra-se no Github [1]

# u_n, v_n -> concentracao dos reagentes no tempo n (matriz Nx x Ny)
# u_n1, v_n1 -> concentracao dos reagentes no tempo n+1

while t < t_max:

    for i in range(Nx):
        i_e = (i - 1) % Nx  # vizinho a esquerda de 0 é o da ultima posicao
        i_d = (i + 1) % Nx  # vizinho a direita da ultima posicao é o zero

        for j in range(Ny):
            j_e = (j - 1) % Ny
            j_d = (j + 1) % Ny

            # calcular u e v no tempo n+1
            u_n1[i, j] = u_n[i, j] + dt * f(u_n[i, j], v_n[i, j], b) \
                         + ku * (u_n[i_e, j] + u_n[i_d, j] + u_n[i, j_e] + u_n[i, j_d] - 4 * u_n[i, j])
            v_n1[i, j] = v_n[i, j] + dt * g(u_n[i, j], v_n[i, j], b) \
                         + kv * (v_n[i_e, j] + v_n[i_d, j] + v_n[i, j_e] + v_n[i, j_d] - 4 * v_n[i, j])

    # atualizar u_n e v_n
    for i in range(Nx):
        for j in range(Ny):
            u_n[i, j] = u_n1[i, j]
            v_n[i, j] = v_n1[i, j]

    t += dt

Resultados

Esta seção se dedica ao resultados das simulações realizadas, que mostram o caráter oscilatório do Brusselator ao longo do tempo e também como se dá a difusão dos reagentes em um recipiente em duas dimensões. As constantes utilizadas foram identificadas abaixo para plotar os gráficos, exceto quando indicado o contrário.

Tabela de Constantes e Valores
Símbolo	Nome	Valor
$N_x$	Dimensão analisada ao longo do eixo x	50
$N_y$	Dimensão analisada ao longo do eixo y	50
$a$	Constante relativa à concentração do reagente A	1
$b$	Constante relativa à concentração do reagente B	1.7
$D_u$	Constante de difusão do reagente U	0.5
$D_v$	Constante de difusão do reagente V	1
$u_0$	Concentração de u no tempo inicial (t=0)	1
$v_0$	Concentração de v no tempo inicial (t=0)	1
$dt$	Passo de tempo entre iterações	0.01
$ds$	Unidade de avanço dos eixos no espaço	1

As posições das condições iniciais foram escolhidas arbitrariamente, tendo um total de 6 modos que serão explicados na sequência, sendo eles:

1. Condição em formato do sinal "+".
2. Condição de borda.
3. Condição aleatória.
4. Condição de nove pontos centrais.
5. Condição da borda completa.

Simulação

Nas figuras abaixo temos a variação das concentrações dos reagentes ao longo do tempo e os diagramas de fase. Para $a=1$ e $b=1.7$ notamos que a solução do sistema converge para $u(t)=1$ e $v(t)=1.7$. Esse resultado está de acordo com o esperado, pois o ponto $(1,1.7)$ é um ponto fixo atrator (satisfaz que $1-b+a^2>0$), conforme discutido na seção de análise de estabilidade do Brusselator. Notamos também, pelo diagrama de fase, que mesmo alterando as condições iniciais do problema, a solução converge para o estado estacionário.

Figura 1: Solução do Brusselator para a = 1 e b = 1.7.

Para $a=1$ e $b=3$, temos que $1-b+a^2<0$, portanto a solução do sistema não converge. O que observamos é um comportamento periódico causado pela existência do ciclo limite.

Figura 2: Solução do Brusselator para a = 1 e b = 3.

Podemos pensar nas funções de concentração $u$ e $v$ como ondas que sofrem interferência destrutiva, pois estão fora de fase. A animação abaixo mostra este comportamento, observe o reagente $U$, representado pela cor amarela, no gráfico da direita. Note que quando o ponto do gráfico da esquerda chega ao mínimo local, o valor de $U$ é alto e o de $V$ muito baixo, fazendo com que o gráfico da direita fique muito amarelo.

Figura 3: Comportamento do reagente.

É interessante pensar no gráfico da esquerda como a análise em um ponto específico dentro do recipiente, por exemplo, se $Nx=Ny=50$ estamos analisando o ponto $(10,10)$. Tendo isso em mente conseguimos analisar os outros 3 cantos da figura, já que, por simetria da posição inicial escolhida, os cantos são iguais.

Figura 4: Análise de um ponto.

---

Condições Iniciais

Realizamos as simulações variando espacialmente as condições iniciais, para verificar os diferentes padrões formados pelo sistema.

Condição em Formato do Sinal "+"

O reagente $u_0$ é distribuído no centro formando um sinal de "+".

Figura 4: Comportamento em formato de "+" para o reagente u. Para o reagente v foi utilizado a condição de borda.

# u_n [Nx, Ny]

# Centro
u_n[int(Nx / 2), int(Ny / 2)] = u0


# Lateral do sinal
#       X                   Y
u_n[int(Nx / 2) + 1, int(Ny / 2)] = u0
u_n[int(Nx / 2) - 1, int(Ny / 2)] = u0


# Altura do sinal
#       X                   Y
u_n[int(Nx / 2), int(Ny / 2) + 1] = u0
u_n[int(Nx / 2), int(Ny / 2) - 1] = u0

---

Condição de borda

O reagente $v_0$ foi distribuído nos 4 cantos do recipiente.

Figura 5: Analisando u com condição de borda "+". Para o reagente v foi utilizado a condição de borda. Foi utilizado b=3 para mostrar a diferença na simulação.

# u_n [Nx, Ny]

# Sup Esq
v_n[0, 0] = v0


# Sup Dir
v_n[0, Nx - 1] = v0


# Inf Esq
v_n[Nx - 1, 0] = v0


# Inf Dir
v_n[Nx - 1, Nx - 1] = v0

---

Condição aleatória

Os reagentes são aleatoriamente distribuídos ao longo do recipiente, a quantidade de focos de reagente também é aleatória para $u_0$ e $v_0$.

Figura 6: Posição aleatória.

from random import *

# u_n [Nx, Ny]

aleatorio = randint(int(Nx / 5)


# Aleatório para u
for _ in range(aleatorio, Nx)):
    u_n[randint(0, size), randint(0, size)] = u0


# Aleatório para v
for _ in range(aleatorio, Nx)):
    v_n[randint(0, size), randint(0, size)] = v0

---

Condição de nove pontos centrais

No centro divide 9 pontos igualmente espaçados de tamanho $N_x/4$, para $u_0$. Foi utilizado para $u_0$ na figura 3.

Figura 7: Utilizando a condição de borda para v.

# Nove pontos centrais (u0)
mid = int(Nx / 2); mov = int(Nx / 4)

# Meio & direita
u_n[mid, mid] = u0
u_n[mid + mov, mid + mov] = u0
# Meio & esquerda
u_n[mid + mov, mid] = u0
u_n[mid, mid + mov] = u0
# Meio, p/baixo & p/lados
u_n[mid + mov, mid - mov] = u0
u_n[mid - mov, mid - mov] = u0
# Meio, p/cima & p/lados
u_n[mid - mov, mid] = u0
u_n[mid, mid - mov] = u0
u_n[mid - mov, mid + mov] = u0

---

Condição da borda completa

Completa toda a borda do recipiente com reagente $v_0$. Foi utilizado para $v_0$ na figura 3.

Figura 8: Utilizando a condição de borda para v e b=3.

# v_n = [Nx, Ny]


# Toda a borda (v0)


for i in range(Nx):

    # Completa primeira coluna e primeira linha
    v_n[0, i] = v0
    v_n[i, 0] = v0

    # Completa última coluna e última linha
    v_n[Nx - 1, i] = v0
    v_n[i, Nx - 1] = v0