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 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. ^[1]. 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^[2]. 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^[3] é 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^[4].

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

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

{\frac {\partial v}{\partial t}}=bu-u^{2}v+D_{v}\left({\frac {\partial ^{2}v}{\partial x^{2}}}+{\frac {\partial ^{2}v}{\partial y^{2}}}\right)

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^{*})=\left(a,{\frac {b}{a}}\right)$ . A matriz jacobiana $J^{*}$ no ponto crítico $(u^{*},v^{*})$ é dada por

J^{*}={\begin{bmatrix}{\frac {\partial f}{\partial u}}{\Big |}_{u^{*},v^{*}}&{\frac {\partial f}{\partial v}}{\Big |}_{u^{*},v^{*}}\\{\frac {\partial g}{\partial u}}{\Big |}_{u^{*},v^{*}}&{\frac {\partial g}{\partial v}}{\Big |}_{u^{*},v^{*}}\end{bmatrix}}={\begin{bmatrix}b-1&a^{2}\\-b&-a^{2}\end{bmatrix}}

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 $\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^[6], 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^[6], foi feita a análise de ponto fixo. O único ponto fixo encontrado é $(u^{*},v^{*})=\left(a,{\frac {b}{a}}\right)$ . Manipulando a matriz Jacobiana nesse ponto e analisando os autovalores, obtém-se que a solução é atraída para o ponto $\left(a,{\frac {b}{a}}\right)$ 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+\Delta t)-f(x,y,t)}{\Delta t}}+{\mathcal {O}}(\Delta t^{2})

{\frac {\partial ^{2}f(x,y,t)}{\partial x^{2}}}\to {\frac {f(x-\Delta x,y,t)-2f(x,y,t)+f(x+\Delta x,y,t)}{(\Delta x)^{2}}}+{\mathcal {O}}(\Delta x^{3})

{\frac {\partial ^{2}f(x,y,t)}{\partial y^{2}}}\to {\frac {f(x,y-\Delta y,t)-2f(x,y,t)+f(x,y+\Delta y,t)}{(\Delta y)^{2}}}+{\mathcal {O}}(\Delta y^{3})

Substituindo nas equações do Brusselator

{\frac {u(x,y,t+\Delta t)-u(x,y,t)}{\Delta t}}=f(u,v)+D_{u}\left({\frac {u(x-\Delta x,y,t)-2u(x,y,t)+u(x+\Delta x,y,t)}{(\Delta x)^{2}}}+{\frac {u(x,y-\Delta y,t)-2u(x,y,t)+u(x,y+\Delta y,t)}{(\Delta y)^{2}}}\right)

{\frac {v(x,y,t+\Delta t)-v(x,y,t)}{\Delta t}}=g(u,v)+D_{v}\left({\frac {v(x-\Delta x,y,t)-2v(x,y,t)+v(x+\Delta x,y,t)}{(\Delta x)^{2}}}+{\frac {v(x,y-\Delta y,t)-2v(x,y,t)+v(x,y+\Delta y,t)}{(\Delta y)^{2}}}\right)

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\Delta t,2\Delta t,3\Delta t,...,t_{max}

x=0,1\Delta x,2\Delta x,3\Delta x,...,N_{x}

y=0,1\Delta y,2\Delta y,3\Delta y,...,N_{y}

Utilizando a notação $f(i\Delta x,j\Delta y,n\Delta t)=f_{i,j}^{n}$ , assumindo $\Delta x=\Delta y=\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})\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})\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}\Delta t}{(\Delta s)^{2}}}$ e $K_{v}={\frac {D_{v}\Delta t}{(\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}\left(\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}\right)+\Delta tf_{u}\delta u_{i,j}^{n}+\Delta tf_{v}\delta v_{i,j}^{n}

\delta v_{i,j}^{n+1}=\delta v_{i,j}^{n}+K_{v}\left(\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}\right)+\Delta tg_{u}\delta u_{i,j}^{n}+\Delta tg_{v}\delta v_{i,j}^{n}

{\delta u^{n} \choose \delta v^{n}}=A^{n}e^{i(k_{x}i\Delta x+k_{y}j\Delta y)}{\delta u_{0} \choose \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 feito por Scholz^[7], obtemos

A_{\pm }={\frac {1}{2}}\left[\beta \pm {\sqrt {\beta ^{2}-4(\Delta t)^{2}a^{2}b}}\right]

onde

\beta =2-4\gamma (K_{u}+K_{v})+\Delta t(b-1-a^{2})

e

\gamma =sin^{2}\left({\frac {k_{x}\Delta s}{2}}\right)+sin^{2}\left({\frac {k_{y}\Delta s}{2}}\right)

Os valores $\Delta t=0.01$ e $\Delta s=1$ satisfazem a condição de estabilidade, como mostrado em Scholz^[7].

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:

Condição em formato do sinal "+".
Condição de borda.
Condição aleatória.
Condição de nove pontos centrais.
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 $\left(1,1.7\right)$ é 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

Referências

↑ J. Tyson, Some further studies of nonlinear oscillations in chemical systems, J. Chem. Phys. 58 (1973) 3919
↑ Turing A. M., The chemical basis of morphogenesis. 1953. Bull. Math. Biol. 52, 153, discussion 119 (1990).
↑ I. Prigogine, R. Lefever, Symmetries breaking instabilities in dissipative systems II. J. Phys. Chem. 48, 1695–1700 (1968)
↑ R. Lefever and G. Nicolis, Chemical instabilities and sustained oscillations, J. Theor. Biol. 30 (1971)
↑ G. Adomian, The diffusion-Brusselator equation. Comput. Math. Appl. 29, 1–3 (1995)
↑ ^6,0 ^6,1 Twizell, E. H., Gumel, A. B., & Cao, Q. (1999). Journal of Mathematical Chemistry, 26(4), 297–316
↑ ^7,0 ^7,1 Scholz, Christian, Morphology of Experimental and Simulated Turing Patterns, 2009, p. 27-29.

[1] J. Tyson, Some further studies of nonlinear oscillations in chemical systems, J. Chem. Phys. 58 (1973) 3919

[2] Turing A. M., The chemical basis of morphogenesis. 1953. Bull. Math. Biol. 52, 153, discussion 119 (1990).

[3] I. Prigogine, R. Lefever, Symmetries breaking instabilities in dissipative systems II. J. Phys. Chem. 48, 1695–1700 (1968)

[4] R. Lefever and G. Nicolis, Chemical instabilities and sustained oscillations, J. Theor. Biol. 30 (1971)

[5] G. Adomian, The diffusion-Brusselator equation. Comput. Math. Appl. 29, 1–3 (1995)

[t-6] 6,0 ^6,1 Twizell, E. H., Gumel, A. B., & Cao, Q. (1999). Journal of Mathematical Chemistry, 26(4), 297–316

[s-7] 7,0 ^7,1 Scholz, Christian, Morphology of Experimental and Simulated Turing Patterns, 2009, p. 27-29.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

Modelo Brusselator de Reação-Difusão

Índice

Modelo de Brusselator

Análise da estabilidade do sistema

Análise de ponto crítico

Análise de ponto fixo

Método FTCS

Análise de estabilidade

Implementação

Resultados

Simulação

Condições Iniciais

Condição em Formato do Sinal "+"

Condição de borda

Condição aleatória

Condição de nove pontos centrais

Condição da borda completa

Referências

Menu de navegação

Modelo Brusselator de Reação-Difusão

Modelo de Brusselator

Análise da estabilidade do sistema

Análise de ponto crítico

Análise de ponto fixo

Método FTCS

Análise de estabilidade

Implementação

Resultados

Simulação

Condições Iniciais

Condição em Formato do Sinal "+"

Condição de borda

Condição aleatória

Condição de nove pontos centrais

Condição da borda completa

Referências

Menu de navegação

Pesquisa