Modelo de Turing

EM CONSTRUÇÃO

Equações de Turing

Simulações computacionais que envolvem equações diferenciais parciais (EDP's) são usualmente modeladas através da discretização das variáveis espaciais e temporais. Algumas dessas equações descrevem comportamentos difusivos no sistema, sendo chamadas de equações de difusão. Tais equações envolvem variáveis de estado que apresentam variações temporal e espacial e coeficientes de difusão no sistema, além de outros parâmetros que influenciam na evolução dos estados. Dentro desse ramo de equações, encontra-se o Modelo de Turing, desenvolvido por Alan Turing, que utiliza como base a concentração de espécies em um sistema, avaliando sua reação, difusão e variação espacial e temporal. São muitas as aplicações do modelo, principalmente em ramos como biologia e química, envolvendo problemas com formação de padrões. A seguir, descrevemos sua formulação matemática.

Sejam $u$ e $v$ as concentrações das espécies que serão analisadas. Sejam $a,b,c$ e $d$ parâmetros e $h$ e $k$ constantes. Os coeficientes de difusão são $D_u$ e $D_v$, cada um associado a uma das concentrações. O Modelo de Turing é dado pelas EDP's

$\frac{\partial u}{\partial t}=a(u-h)+b(v-k)+D_u{\mathrm{\nabla }}^{2}u$

$\frac{\partial v}{\partial t}=c(u-h)+d(v-k)+D_v{\mathrm{\nabla }}^{2}v$

Note que certa parte de cada equação depende apenas dos parâmetros e das concentrações. Podemos, portanto, utilizar funções de variáveis $u$ e $v$ para descrever o sistema, de modo que

$u_t=D_u{\mathrm{\nabla }}^{2}u+f(u,v)$

$v_t=D_v{\mathrm{\nabla }}^{2}v+g(u,v)$

Estabilidade e Instabilidade no Modelo de Turing

Pontos de Equilíbrio

Vimos que o modelo de Turing depende de parâmetros $a,b,c,d$, de constantes $h$ e $k$ e dos coeficientes de difusão.

Afirmação: Se $D_u=D_v=0$, temos ($v_{eq},u_{eq})=(h,k)$ como o único ponto de equilíbrio.

Demonstração: Mostraremos que $(h,k)$ é ponto de equilíbrio. De fato, ao aplicarmos esse ponto na equação do modelo de Turing, temos

$u_t{|}_{u=k}=0=v_t{|}_{v=h}$

para mostrar que é único, suponha que existem dois pontos de equilíbrio, a saber, $(v_1,u_1)$ e $(v_2,u_2)$. Vemos que, como as equações diferenciais em cada ponto fixo são iguais a zero, temos

$a(u_1-u_2)+b(v_1-v_2)=0$

$c(u_1-u_2)+d(v_1-v_2)=0$

Consequentemente, devemos ter

$(b-\frac{ad}{c})(v_1-v_2)=0⟹v_1=v_2$.

Do mesmo modo, $u_1=u_2$. Portanto, o ponto de equilíbrio é único nessas circunstâncias.

Estabilidade de Sistemas Reativos-Difusivos

Para estudarmos a estabilidade dos sistemas reativos-difusivos precisamos encontrar os autovalores da matriz

$(J-D{\omega }^2){|}_{f=f_{eq}}$

Onde $J$ é a matriz jacobiana dos termos de reação, $D$ é a matriz diagonal dos termos de difusão e $\omega$ é o parâmetro que determina a frequência espacial das perturbações. A prova dessa afirmação pode ser encontrada nas referências. Aplicando isso ao modelo de Turing obtemos

$(\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)-\left(\begin{array}{cc}{D}_u& 0\\ 0& D_v\end{array}\right){\omega }^2){|}_{(u,v)=(h,k)}=\left(\begin{array}{cc}a-D_u{\omega }^2& b\\ c& d-D_v{\omega }^2\end{array}\right)$

Para esta matriz ser estável precisamos que o determinante dessa matriz seja positivo e o traço seja negativo. Obtemos então

$0<(a-D_u{\omega }^2)(d-D_v{\omega }^2)-bc$

$0>a-D_u{\omega }^2+d-D_v{\omega }^2$

Podemos ver que em ambas desigualdades aparecerá o Determinante e o Traço da matriz dos coeficientes de reação. Reescrevendo então obtemos

$a{D}_v{\omega }^2+d{D}_u{\omega }^2-D_u{D}_v{\omega }^4<\text{det}(J)$

$D_u{\omega }^2+D_v{\omega }^2>\text{Tr}(J)$

Se o sistema fosse originalmente estável, isto é, $\text{Tr}(J)<0$ e $\text{det}(J)>0$ a segunda desigualdade é sempre verdade, mas a primeira não. Podemos ver que a desigualdade pode ser violada se

$g(z)=-D_u{D}_v{z}^2+(a{D}_v+d{D}_u)z-\text{det}(J)$

tomar valores positivos para algum $z>0$ (Onde $z={\omega }^2$). Podemos reescrever $g(z)$ da forma

$g(z)=-D_u{D}_v{(z-\frac{a{D}_v+d{D}_u}{2D_u{D}_v})}^2+\left(\frac{(a{D}_v+d{D}_u{)}^2}{4D_u{D}_v}\right)-\text{det}(J)$

Vamos então analisar duas opções onde podemos ter $z>0$: (1) o ponto mais alto da função fica no lado positivo do eixo $z$ ou (2) o ponto mais alto da função fica no no lado negativo do eixo $z$. Na primeira opção (1) a única condição é que o máximo da função seja acima do eixo $z$, logo teremos

$\frac{(a{D}_v+d{D}_u{)}^2}{4D_u{D}_v}-\text{det}(J)>0$

Se o máximo da função estiver no lado negativo de $z$, opção (2), a condição é que o ponto que interceptar $g(z)$ deve ser positivo, podemos escrever isso da forma

$g(0)=-\text{det}(J)>0$

Podemos ver que a opção (2) nunca pode ser obtida, pois o modelo é inicialmente estável e $\text{det}(J)>0$, assim a única opção para que a difusão desestabilize o sistema é quando

$a{D}_v+d{D}_u>2\sqrt{D_u{D}_v\text{det}(J)}$

e com essa desigualdade é fácil checar quais valores dos parâmetros $a$, $b$, $c$ e $d$ e dos coeficientes de difusão $D_u$ e $D_v$ formarão padrões de Turing.

Implementação

Para resolver numericamente as equações de Turing iremos utilizar o método FTCS (Forward Time Central Space). O método FTCS é o mais simples e consiste em discretizar a derivada em $t$ de forma não simetrizada. Obtemos as seguintes discretizações para uma função genérica $f(r,t)$

$\frac{\partial f}{\partial t}=\frac{f(r,t+dt)-f(r,t)}{dt}$

$\frac{{\partial }^{2}f}{\partial r^2}=\frac{f(r+dr,t)-2f(r,t)+f(r-dr,t)}{d{r}^2}$

Onde $\overrightarrow{r}$ é o vetor posição, que neste trabalho utilizamos apenas duas dimensões, $\overrightarrow{r}=(x,y)$.

Podemos discretizar as equações de Turing diretamente com o método FTCS. Talvez o único problema seja o laplaciano, porém basta escrever da forma

${\mathrm{\nabla }}^{2}f(r,t)=\frac{{\partial }^{2}f}{\partial x^2}+\frac{{\partial }^{2}f}{\partial y^2}$

Assim podemos utilizar a discretização simetrizada e obter

${\mathrm{\nabla }}^{2}f(r,t)=\frac{f(x+dx,t)-2f(r,t)+f(x-dx,t)}{d{x}^2}+\frac{f(y+dy,t)-2f(r,t)+f(y-dy,t)}{d{y}^2}$

Ao tomarmos $dx=dy=dh$, que faremos aqui, podemos simplificar a discretização do laplaciano para

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

Então obtemos que as equações de Turing discretizadas pelo método FTCS, em notação discreta, são dadas por

$u_{i,j}^{n+1}=u_{i,j}^n+[a(u_{i,j}^n-h)+b(v_{i,j}^n-k)+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}{d{h}^2}]dt$

$v_{i,j}^{n+1}=v_{i,j}^n+[c(u_{i,j}^n-h)+d(v_{i,j}^n-k)+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}{d{h}^2}]dt$

Onde $i$ e $j$ são os índices espaciais e $n$ é o índice temporal.

Utilizamos uma rede quadrada de tamanho $L\times L$ com condições de contorno periódicas. O sistema inicia próximo do equilibrio $(u,v)=(h,k)$ e então é aplicado um pequeno ruído para começar a difusão. O ruído é muito importante, sem ele o sistema ficaria sempre no equilíbrio. O ruído também deve ser pequeno suficiente para quebrar o estado inicial, mas não grande suficiente para causar instabilidades numéricas na simulação. O ruído utilizado aqui consiste em números aleatórios no intervalo $[-0,003:0,003]$. Tomamos $dh=1/L$.

Resultados e Discussão

Na rede $100\times 100$ escolhendo $(a,b,c,d)=(1,-1,2,-1.5)$ e $(D_u,D_v)=(10^{-4},6\times 10^{-4})$, podemos verificar se a desigualdade é respeitada.

$a{D}_v+d{D}_u=4.5\times 10^{-4}$

$2\sqrt{D_u{D}_v\text{det}(J)}\approx 3.464\times 10^{-4}$

A desigualdade é sim respeitada. Abaixo podemos ver a simulação das concentrações de $u$ e de $v$.

Modelo Turing com $(a,b,c,d)=(1,-1,2,-1.5)$ e $(D_u,D_v)=(10^{-4},6\times 10^{-4})$
Simulação do Modelo de Turing da concentração de $u$.	Simulação do Modelo de Turing da concentração de $v$.

Quanto mais amarelo maior a densidade, assim como quanto mais roxo menor a densidade. Outra escolha de constantes é $(a,b,c,d)=(0.5,-1,0.75,-1;)$ e $(D_u,D_v)=(10^{-4},10^{-3})$. Verificando que a desigualdade é respeitada:

$a{D}_v+d{D}_u=4\times 10^{-4}$

$2\sqrt{D_u{D}_v\text{det}(J)}\approx 3.162\times 10^{-4}$

Novamente respeitamos a desigualdade. Na rede $100\times 100$ obtemos a seguinte simulação das concentrações de $u$ e de $v$.

Modelo de Turing com $(a,b,c,d)=(0.5,-1,0.75,-1;)$ e $(D_u,D_v)=(10^{-4},10^{-3})$
Simulação do Modelo de Turing da concentração de $u$.	Simulação do Modelo de Turing da concentração de $v$.

É possível ver que o sistema rapidamente forma um padrão, onde existem algum aglomerados com maior densidade e diversos pontos onde existe baixa densidade. Diferentes conjuntos das constantes e dos coeficientes de difusão formam diferentes padrões.

Podemos estudar também a estudar a razão crítica dos coeficientes de difusão, onde o sistema passa a formar padrões. Definindo um novo parâmetro $\rho =D_v/D_u$ podemos rearranjar a desigualdade necessária para formar padrões da forma

$a\rho D_u+d{D}_v>2\sqrt{\rho D_u^2\text{det}(J)}$

$a\rho +d>2\sqrt{\rho \text{det}(J)}$

Por exemplo, aplicando isso para um conjunto de constantes que simulamos $(a,b,c,d)=(1,-1,2,-1.5)$ chegaremos a

$\rho -1.5>2\sqrt{0.5\rho }$

$\rho >4.5$

Isso significa que a difusão de $v$ deve ser quase cinco vezes maior que a difusão de $u$ para gerar os padrões. A difusão de $u$ age mais localmente enquanto a de $v$ age em áreas maiores. Uma transição de fase é observada perto de $\rho =4.5$. Essa transição de fase entre homogenização e formação de padrões é chamada de Bifurcação de Turing .

Programas Utilizados

Programas na linguagem C. Para utilizar os programas, abra o terminal e compile da forma

$ gcc prog.c -lm

Onde prog.c é o programa que deseja utilizar. E execute da seguinte maneira

$ ./a.out PARAMETROS

onde o segundo termo são os parâmetros do sistema, argumento dos programas. Dentro dos programas tem um exemplo de execução. O programa da simulação possuei uma diretiva de compilação para visualização do sistema ao decorrer da execução. Para utilizar é necessário ter o gnuplot instalado e compilar da forma

$ gcc -DGNU prog.c -lm

e então executar da maneira

$ ./a.out PARAMETROS | gnuplot

Estabilidade das Constantes

Simulação do Modelo de Turing

Conclusão

Partindo das equações diferenciais que caracterizam o modelo de Turing, analisamos os pontos de estabilidade no caso específico $D_u=D_v=0$ e no caso mais geral, utilizando os autovalores da matriz $(J-D{\omega }^2){|}_{f=f_{eq}}$. Com isso, foi obtido uma equação que determina se os parâmetros desejados formarão padrões de Turing. Descrevemos a implementação de um algoritmo que aplica o método FTCS para encontrar soluções numericamente. Para alguns parâmetros selecionados, foram feitas simulações do modelo para a concentração de $u$ e $v$, nas quais notou-se a formação de padrões que dependem dos parâmetros. Acima, foram deixados os programas computacionais utilizados na resolução numérica do método de Turing.