Modelo Brusselator de Reação-Difusão

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

Figura 1: Variação da concentração em um recipiente de 100x100. Áreas em amarelo correspondem a uma grande quantidade de reagente U.

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 1970 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. [J. Tyson, Some further studies of nonlinear oscillations in chemical systems, J. Chem. Phys. 58 (1973) 3919.] Turing [ref?] 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 (1970) é dado por [I. Prigogine, R. Lefever, Symmetries breaking instabilities in dissipative systems II. J. Phys. Chem. 48, 1695–1700 (1968)]:

$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 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 “limit-cycle oscilations”. [ R. Lefever and G. Nicolis, Chemical instabilities and sustained oscillations, J. Theor. Biol. 30 (1971) 267.].

As equações diferenciais parciais associadas com o sistema Brusselator são dadas por(G. Adomian, The diffusion-Brusselator equation. Comput. Math. Appl. 29, 1–3 (1995)):

$\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$ constantes de difusão.

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 $\partial x=\partial y=0$:

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

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

Onde $u=u(t)$ e $v=v(t)$ e $a$ e $b$ são constantes positivas e reais. A matriz jacobiana $J^{*}$ no ponto crítico $(u^{*},v^{*})$ é dada por

$J^{*}=\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 ou determinam a existência de um ciclo limite. As propriedades de estabilidade ou a existência de um ciclo limite estão sumarizadas na tabela abaixo, em relação a figura 1.

[[Arquivo:]]

Utilizando a teoria Hopf, é mostrado que o ponto crítico perde sua estabilidade quando A e B movem da região 2 para região 3, na figura 1, atravessando a curva $1-b+a^2=0$. Uma bifurcação Hopf ocorre quando ao passo que essa curva é atravessada e um ciclo limite estável existe para A e B nas regiões 1 e 2, mas não para A e B nas regiões 3 e 4.

Conclui-se que a curva $1-b+a^2=0$ governa a estabilidade do sistema.

Análise de ponto fixo

O estado estacionário do sistema pode ser encontrado igualando o coeficiente de difusão, e portanto as derivadas parciais, a zero. Percebe-se que esse sistema converge para os pontos fixos:

$u=a$

$v=\frac{b}{a}$

Em [twizell] a manipulação da matriz jacobiana nos pontos fixos resulta nos seguintes autovalores

${\mu }_1=\frac{1+\frac{1}{2}\ell (1-b+a^2)-\frac{1}{4}{\ell }^2{a}^2}{1+\frac{1}{2}\ell (1-b+a^2)+\frac{1}{4}{\ell }^2{a}^2}$

${\mu }_2=\frac{1-\frac{1}{2}\ell (1-b+a^2)-\frac{1}{4}{\ell }^2{a}^2}{1+\frac{1}{2}\ell (1-b+a^2)+\frac{1}{4}{\ell }^2{a}^2}$

Onde fica claro que os denominadores dos autovalores são sempre positivos quando $1-b+a^2>0$ and $\ell >0$. As inequações

$|{\mu }_1|<1$ e $|{\mu }_2|<1$

São verdadeiras sempre que $1-b+a^2=0$. Portanto, uma condição suficiente para o ponto fixo $(a,\frac{b}{a})$ atrair a sequência gerada pelo sistema é $1-b+a^2>0$.

Ainda em [Twizell] o modelo reação-difusão Brusselator foi discretizado e a análise dos pontos fixos concluiu que o sistema converge para $u(x,t)=a$ e $v(x,t)=\frac{b}{a}$, sendo esse o único estado estacionário do sistema.

Concluiu que também que o sistema apresenta estado oscilatório quando

$1-b+a<0$

Estado em que o sistema não converge para nenhum ponto.

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 do método

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})=A^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 $|A|<1$. Após inserir as soluções nas equações linearizadas e realizar manipulações algébricas, como em [ref] obtemos

$A_{\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 [ref].

Resultados

Esta seção se dedica ao resultados da simulação, utilizando valores de constantes de outros artigos como o [Nome do artigo e referência]. As constantes utilizadas foram identificadas na próxima subseção, com exceção nos gráficos em que indicamos outros valores. As simulações mostram o caráter oscilatório do Brusselator ao longo do tempo e também como se da a reação e difusão em um recipiente em duas dimensões. O código utilizado (incompleto) estará logo após a seção de resultados.

Definindo Constantes

Durante as simulações em 2D foram utilizados no código:

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.1
$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).	2
$dt$	Passo de tempo entre iterações.	0.1
$ds$	Unidade de avanço dos eixos no espaço.	1

Constantes que dependem diretamente de outros parâmetros não recebem um valor fixo, como é o caso de $k_u$ e $k_v$ que são constantes de estabilidade. Utilizamos estes valores pois foi utilizado no artigos. Já as posições das condições iniciais foram escolhidas arbitrariamente, tendo um total de 6 modos que serão explicados logo à frente, 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

No gráfico abaixo podemos notar o caráter ondulatório da variação da concentração ao longo do tempo, as linhas (azul e vermelho) significam a taxa de concentração no recipiente.

Note que quando $U$ diminui, $V$ aumenta rapidamente, mantendo-se um ciclo oscilatório em que a amplitude da onda de $U$ e $V$ vão diminuindo ao longo do tempo, como mostrado nas figuras.

Figura 1: Gráfico de concentrações.

Pensando no gráfico da reação-difusão como ondas, é correto afirmar que ambas juntas formam um padrão de interferência destrutivo. A animação abaixo mostra exatamente este comportamento, observando o reagente $U$ 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 (amarelo significa uma alta concentração de $U$ sendo formada).

Figura 2:: Comportamento do reagente.

---

Condições Iniciais

Para observar melhor os resultados, foram feitas diversas condições iniciais que foram separadas em tópicos citados em Definindo Constantes.

Levando em consideração que u_n é um array de tamanho $N_x$ por $N_y$.

Condição em Formato do Sinal "+"

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

# 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.

# 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$.

# Aleatório para u
for _ in range(randint(int(Nx / 5), Nx)):
    u_n[randint(0, size), randint(0, size)] = u0
# Aleatório para v
for _ in range(randint(int(Nx / 5), Nx)):
    v_n[randint(0, size), randint(0, size)] = v0

Condição de nove pontos centrais

Foi utilizado na figura 2.

No centro divide 9 pontos igualmente espaçados de tamanho $N_x/4$, para $u_0$.

# Nove pontos centrais (u0)
mid = int(Nx / 2)  # Centro
mov = int(Nx / 4)  # Movendo para cima e para os lados

u_n[mid, mid] = u_n[mid + mov, mid + mov] = u_n[mid + mov, mid] = u_n[mid, mid + mov] = u_n[mid + mov, mid - mov] = u0
u_n[mid - mov, mid - mov] = u_n[mid - mov, mid] = u_n[mid, mid - mov] = u_n[mid - mov, mid + mov] = u0

Condição da borda completa

Foi utilizado na figura 2.

Completa toda a borda do recipiente com reagente $v_0$.

# Toda a borda (v0)
for i in range(Nx):
    v_n[0, i] = v0
    v_n[i, 0] = v0
    v_n[Nx - 1, i] = v0
    v_n[i, Nx - 1] = v0

Implementação

O método foi implementado em Python, considerando $\mathrm{\Delta }t=0.1,\mathrm{\Delta }s=1,N_x=N_y=50$, variando as constantes $a$ e $b$ e as condições iniciais do problema.

• Código completo no GitHub[1]

import numpy as np

# Constantes
Nx = Ny = 25
a = 1
b = 1.7
Du = 0.1
Dv = 1
t_max = 40
dt = 0.1
ds = 1

ku = Du * dt / (ds ** 2)
kv = Dv * dt / (ds ** 2)

# Valores iniciais
u0 = 1
v0 = 2
t = 0


def f(u, v):
    return a - (b + 1) * u + u * u * v



def g(u, v):
    return b * u - u * u * v


# vetores no tempo n
u_n = np.zeros((Nx, Ny))
v_n = np.zeros((Nx, Ny))

# vetores no tempo n+1
u_n1 = np.zeros((Nx, Ny))
v_n1 = np.zeros((Nx, Ny))

# Condições iniciais (u e v distribuidos aleatoriamente no espaço
for _ in range(randint(int(Nx / 5), Nx)):
    u_n[randint(0, size), randint(0, size)] = u0
for _ in range(randint(int(Nx / 5), Nx)):
    v_n[randint(0, size), randint(0, size)] = v0

# FTCS
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

            u_n1[i, j] = u_n[i, j] + dt * f(u_n[i, j], v_n[i, j]) + 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]) + 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