Exibir código-fonte para 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. <ref>J. Tyson, Some further studies of nonlinear oscillations in chemical systems, J. Chem. Phys. 58 (1973) 3919</ref>. 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<ref>Turing A. M., The chemical basis of morphogenesis. 1953. Bull. Math. Biol. 52, 153, discussion 119 (1990).</ref>. 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<ref>I. Prigogine, R. Lefever, Symmetries breaking instabilities in dissipative systems II. J. Phys. Chem. 48, 1695–1700 (1968)</ref> é dado por: 

:<math>A \to U</math>         (1.a)
:<math>B + U \to V + D</math> (1. b)
:<math>2U + V \to 3U</math> (1.c)
:<math>U \to E</math> (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<ref>R. Lefever and G. Nicolis, Chemical instabilities and sustained oscillations, J. Theor. Biol. 30 (1971)</ref>.

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

:<math>\frac{\partial u}{\partial t} = a + u^2v - (b+1)u + D_u \left(\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2} \right)</math>
:<math>\frac{\partial v}{\partial t} = bu - u^2v + D_v \left(\frac{\partial ^2 v}{\partial x^2} + \frac{\partial ^2 v}{\partial y^2}\right)</math>

onde <math>u(x, y, t)</math> e <math>v(x, y, t)</math> são as concentrações a serem investigadas em função de tempo e espaço, <math>a</math> e <math>b</math> são constantes relativas às concentrações dos reagentes A e B, e <math>D_u</math> e <math>D_v</math> 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 <math>D_u = D_v = 0</math>:

:<math> \frac{du}{dt} = f(u, v) = a + u^2v - (b+1)u </math>

:<math> \frac{dv}{dt} = g(u, v) = bu - u^2v </math>

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

:<math>
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}
</math>


Os autovalores de <math>J^*</math> são os valores <math>\lambda</math> que satisfazem a equação caracterísitca

:<math> \lambda ^2 + (1 - b + a^2)\lambda + a^2 = 0 </math>

Os autovalores claramente mostram dependência em <math>1-b+a^2</math> e no determinante <math>\Delta = (1 - b + a^2)^2 - 4a^2</math>. Esses autovalores governam a estabilidade do ponto crítico e determinam a existência de um ciclo limite. Conforme discutido em Twizell<ref>Twizell, E. H., Gumel, A. B., & Cao, Q. (1999). Journal of Mathematical Chemistry, 26(4), 297–316 </ref>, quando <math>Re(\lambda) < 0 </math> o sistema é estável e não existe ciclo limite. Já para <math>Re(\lambda) > 0 </math> o sistema é instável e existe ciclo limite.

===Análise de ponto fixo===

Ainda em Twizell<ref>Twizell, E. H., Gumel, A. B., & Cao, Q. (1999). Journal of Mathematical Chemistry, 26(4), 297–316 </ref>, foi feita a análise de ponto fixo. O único ponto fixo encontrado é <math>(u^*, v^*) = \left ( a, \frac{b}{a} \right)</math>. Manipulando a matriz Jacobiana nesse ponto e analisando os autovalores, obtém-se que a solução é atraída para o ponto <math>\left ( a, \frac{b}{a} \right)</math> se <math>1-b+a^2 > 0</math>. 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:

:<math>
\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)
</math>

:<math>
\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)
</math>

:<math>
\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)
</math>


Substituindo nas equações do Brusselator

:<math>
\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)
</math>

:<math>
\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)
</math>

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


Utilizamos discretização do tipo
:<math>t = 0, 1\Delta t, 2\Delta t, 3\Delta t, ..., t_{max}</math>
:<math>x = 0, 1\Delta x, 2\Delta x, 3\Delta x, ..., N_x</math>
:<math>y = 0, 1\Delta y, 2\Delta y, 3\Delta y, ..., N_y</math>


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

:<math>
    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 - 4 u_{i, j}^n)
</math>

:<math>
    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 - 4 v_{i, j}^n)
</math>


onde <math>K_u = \frac{D_u \Delta t}{(\Delta s)^2}</math> e <math>K_v = \frac{D_v \Delta t}{(\Delta s)^2}</math>.

===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: 

:<math>
\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 t f_u \delta u_{i, j}^n + \Delta t f_v \delta v_{i, j}^n
</math>

:<math>
\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 t g_u \delta u_{i, j}^n + \Delta t g_v \delta v_{i, j}^n
</math>

:<math>
{\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}
</math>

Para que o método seja estável, é preciso que <math>|A| < 1</math>. Após inserir as soluções nas equações linearizadas e realizar manipulações algébricas, como feito por Scholz<ref>Scholz, Christian, Morphology of Experimental and Simulated Turing Patterns, 2009, p. 27-29.</ref>, obtemos
:<math>
A_\pm = \frac{1}{2} \left [ \beta \pm \sqrt{\beta ^2 - 4 (\Delta t)^2 a^2 b} \right ]
</math>
onde
:<math>\beta = 2 - 4 \gamma (K_u + K_v) + \Delta t (b - 1 - a^2)</math>
e
:<math>\gamma = sin^2\left( \frac{k_x \Delta s}{2} \right) + sin^2\left( \frac{k_y \Delta s}{2} \right)</math>

Os valores <math>\Delta t = 0.01</math> e <math>\Delta s = 1</math> satisfazem a condição de estabilidade, como mostrado em Scholz<ref>Scholz, Christian, Morphology of Experimental and Simulated Turing Patterns, 2009, p. 27-29.</ref>.

===Implementação===

O método foi implementado em Python, utilizando as constantes definidas na tabela abaixo. 

Os valores de <math>b</math> e as condições iniciais foram variados, como descrito na seção de resultados.

A seguir, o trecho do código relativo à implementação do método. O código completo encontra-se no Github [https://github.com/RicGary/MetCompC/tree/main/Brusselator]

<source lang="python">
# 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
# Nx, Ny -> Tamanho do recipiente

    while t < t_max:

        for i in range(Nx):
            i_e = (i - 1) % Nx  # garante que o vizinho a esquerda de 0 é o da ultima posicao
            i_d = (i + 1) % Nx  # garante que 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], 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
 
</source>

==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).

{| class="wikitable" style="margin-left: auto; margin-right: 0px;" style="text-align:center;" 
|- style="text-align:left;"
! colspan="6" style="text-align:center; font-weight:bold; font-size:medium;" | Tabela de Constantes e Valores
|- style="font-weight:bold;"
| Símbolo
| Nome
| Valor
|-
| <math>N_x</math>
| style="text-align:left;" | Dimensão analisada ao longo do eixo x
| 50
|-
| <math>N_y</math>
| style="text-align:left;" | Dimensão analisada ao longo do eixo y
| 50
|-
| <math>a</math>
| style="text-align:left;" | Constante relativa à concentração do reagente A
| 1
|-
| <math>b</math>
| style="text-align:left;" | Constante relativa à concentração do reagente B
| -
|-
| <math>D_u</math>
| style="text-align:left;" | Constante de difusão do reagente U
| 0.5
|-
| <math>D_v</math>
| style="text-align:left;" | Constante de difusão do reagente V
| 1
|-
| <math>u_0</math>
| style="text-align:left;" | Concentração de u no tempo inicial (t=0)
| -
|-
| <math>v_0</math>
| style="text-align:left;" | Concentração de v no tempo inicial (t=0)
| -
|-
| <math>dt</math>
| style="text-align:left;" | Passo de tempo entre iterações
| 0.01
|-
| <math>ds</math>
| style="text-align:left;" | 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 logo à frente, 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 <math>a = 1</math> e <math>b = 1.7</math> notamos que a solução do sistema converge para <math>u(t) = 1</math> e <math>v(t) = 1.7</math>. Esse resultado está de acordo com o esperado, pois o ponto <math>\left(1, 1.7\right)</math> é um ponto fixo atrator (satisfaz que <math>1-b+a^2 > 0</math>), conforme demonstrado 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.


[[Arquivo:Sol_bru_b1.png|thumb|center|790px|'''Figura 1''': Gráfico de concentrações.]]


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


[[Arquivo:Sol_bru_b3.png|thumb|center|790px|'''Figura 2''': Gráfico de concentrações.]]


Note que quando <math>U</math> diminui, <math>V</math> aumenta rapidamente, mantendo-se um ciclo oscilatório em que a amplitude das ondas de <math>U</math> e <math>V</math> vão diminuindo ao longo do tempo (ou não, dependendo das constantes <math>a</math> e <math>b</math> teremos outros comportamentos), como mostrado nas figuras.

Pensando no gráfico da reação-difusão como ondas, é correto afirmar que ambas juntas formam um padrão de interferência destrutiva. A animação abaixo mostra exatamente este comportamento, observe o reagente <math>U</math> no gráfico da direita. Note que quando o ponto do gráfico da esquerda chega ao mínimo local, o valor de <math>U</math> é alto e o de <math>V</math> muito baixo, fazendo com que o gráfico da direita fique muito amarelo (amarelo significa uma alta concentração de <math>U</math> sendo formada).


[[Arquivo:AnaliseMid.gif|thumb|center|1300px|'''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 <math>Nx=Ny=50</math> estamos analisando o ponto <math>(10,10)</math>. 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.


[[Arquivo:ZoomGraph.jpg|thumb|center|800px|'''Figura 4:''': Análise de um ponto.]]

----

===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 <math>N_x</math> por <math>N_y</math>.

====Condição em Formato do Sinal "+"==== 

O reagente <math>u_0</math> é distribuído no centro formando um sinal de "+".

[[Arquivo:CruzBorda.gif|thumb|right|338px|'''Figura 4:''': Comportamento em formato de "+" para o reagente u. Para o reagente v foi utilizado a condição de borda.]]

<source lang="python">

# 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





</source>

----

====Condição de borda====

O reagente <math>v_0</math> foi distribuído nos 4 cantos do recipiente.

[[Arquivo:CruzBordaB3.gif|thumb|right|338px|'''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.]]

<source lang="python">

# 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

</source>

----

====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 <math>u_0</math> e <math>v_0</math>.

[[Arquivo:AleatorioResize.gif|thumb|right|360px|'''Figura 6:''': Posição aleatória.]]

<source lang="python">

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

</source>

----

====Condição de nove pontos centrais====

Foi utilizado na figura 3.

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

[[Arquivo:9PontosBorda.gif|thumb|right|384px|'''Figura 7:''': Utilizando a condição de borda para v.]]

<source lang="python">

# 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

</source>

----

====Condição da borda completa====
Foi utilizado na figura 3.

Completa toda a borda do recipiente com reagente <math>v_0</math>.

[[Arquivo:9PontosB3.gif|thumb|right|359px|'''Figura 8:''': Utilizando a condição de borda para v e b=3.]]

<source lang="python">

# 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

</source>

==Referências==

<references/>