Modelo de reação-difusão

O Modelo de Reação-Difusão é um modelo matemático teorético de como são formados e regulados padrões característicos durante o desenvolvimento de embriões animais. É conhecido também como Modelo Turing em referência ao cientista e matemático Alan Turing, que iniciou esta linha de pesquisa com seu último artigo publicado antes de sua morte.

A tentativa de modelar fenômenos biológicos matematicamente vêm da dificuldade de se obter relações intuitivas de causa e consequência originários destes fenômenos, frequentemente mediados por múltiplas vias moleculares redundantes ou auto-regulatórias. Assim, estabeleceram-se duas iniciativas frente aos problemas de modelagem:^[1]

1. Analisar quantitativamente todos os componentes de uma rede molecular e simular em computadores suas interações. Esta prática é útil e efetiva em sistemas menos robustos, como vias simples em células únicas e é utilizado extensamente no campo da biologia de sistemas.
2. Omitir informações do sistema propositalmente a fim de simplificá-lo em um modelo que consegue gerar previsões efetivas para um fenômeno complexo mesmo sem contar com todas suas partes. Esta estratégia é requerida quando padrões espaçotemporais são envolvidos na equação, tornando previsões computacionais pouco confiáveis.

Utilizando-se desta segunda estratégia, Alan Turing desenvolveu o Modelo de Reação-Difusão, capaz de descrever como padrões espaciais podem desenvolver-se autonomamente em embriões animais.

Trabalhando sobre o modelo de Reação-Difusão de Turing, cientistas como Crick e Driever e suas respectivas equipes aprimoraram e formalizaram o modelo entre a década de 70 e a década de 90^[2][3]. Recentemente, em 2007, Gregor et. al. aperfeiçoaram o modelo e denominaram-no o Modelo Síntese-Difusão-Degradação (SDD)^[4]. Este modelo consta em um sistema de (ao menos) duas substâncias que difundem-se pelo embrião, regulando sua própria produção e interagindo entre si. A inovação de Turing que permitiu a evolução deste modelo foi a concepção de que uma substância não teria efeito morfógeno antes de interagir com outra substância. Esta quebra de paradigma - a introdução da reação entre moléculas difundidas - permitiu que padrões muito mais complexos do que os ditados pela difusão simples de substâncias em diferentes pontos pré-determinados, dependentes portanto de um pré-padrão genético.

Em sua publicação original, Turing previu que duas moléculas teoréticas poderiam formar 6 estados estáveis:^[5]

1. O sistema converge para um estado uniforme e estável
2. Ocorre uma oscilação uniforme da concentração de morfógenos, como se observa em células musculares cardíacas^[6] ou no ciclo circadiano^[7].
3. O sistema forma um padrão de "ilhas", como por exemplo quando células diferenciadas inibem suas vizinhas de diferenciar. Isto é observado por exemplo nas células diferenciadas neuroprogenitoras do epitélio de embriões de Drosophila melanogaster.^[8]
4. Um estado ainda não observado em sistemas biológicos no qual um sistema descrito em 3 tem padrão oscilante. Requereria ao menos 3 substâncias morfógenas.^[9]
5. Uma onda de gradiente pulsa pelo embrião, como observado na formação de padrões espirais na ameba social Dictyostelium discoideum em agregação^[10] e na onda de cálcio que atravessa o ovo do anuro Xenopus laevis como consequência da entrada de um espermatozóide^[11].
6. Ocorrem padrões estacionários, com comprimento de onda finito. Estes são a grande descoberta de Turing e apropriadamente intitulados de Padrões de Turing. Este padrão de onda não-linear e seu comprimento de onda é regulado por um equilíbrio entre diversas propriedades do sistema, como velocidade de produção das moléculas envolvidas, velocidade de difusão, velocidade de reação e velocidade de degradação destas. A capacidade de auto-regeneração do padrão após distúrbios induzidos em laboratório corroboram a robustez e validade do modelo nas observações.^[12][13]

Para um morfógeno qualquer num embrião de eixo de dimensão L, temos que a função de sua concentração é dada por "c(x,t)", a ver:

${\displaystyle {\frac {\partial c(x,t)}{\partial t}}=s(x,t)+D\bigtriangledown ^{2}c(x,t)-k_{\text{deg}}c(x,t)\mid 0<x<L}$

Para a qual:

s(x,t) representa a função de produção de morfógeno;

D é o coeficiente de difusão;

ᐁ é o vetor gradiente;

k_deg é a taxa de degradação.

Costuma-se utilizar condições de contorno de Neumann^[14] e, portanto, há ausência de fluxo nos extremos. Assim; ${\displaystyle \partial c/\partial x\mid _{\text{x=0,L}}=0}$

Em estabilidade, o gradiente C(x) é delimitado por um equilíbrio entre taxas de difusão e de degradação.^[15][3] O tempo requerido para se alcançar uma situação de equilíbrio no entanto depende da meia-vida do morfógeno e - se houver - mudanças na taxa de síntese.^[16]

Podemos simplificar o sistema acima ao assumirmos algumas idealizações, que são aplicáveis apenas a casos específicos:

1. A síntese ocorre em taxa constante.
2. A síntese ocorre exclusivamente no ápice anterior do embrião.
3. A concentração do morfógeno no ápice posterior do embrião é negligenciável.

Deste modo, as duas primeiras condições impõe que ${\displaystyle \partial c/\partial x\mid _{\text{x=0}}=-p/D}$ , com "p" sendo uma constante da taxa de produção.

Junto com a condição 3, podemos estabelecer o equilíbrio da concentração do morfógeno, que toma a forma de uma característica função de decaimento exponencial ${\displaystyle C(x)\propto e^{-x/\lambda }}$ com comprimento de onda ${\displaystyle \lambda ={\sqrt {D/k_{\text{deg}}}}}$^[15][2].^[3][17].

Para entender matematicamente a instabilidade do sistema proposto, podemos utilizar uma versão simples do modelo, conhecida como Modelo de Reação-Difusão tipo Ativador-Inibidor. neste modelo, um Ativador (A) estimula a produção de um Inibidor (I), que bloqueia a síntese de A e decai com o tempo. As concentrações de A e I respeitam às taxas respectivas:

${\displaystyle {\partial a \over \partial t}=f_{a}A+f_{i}I+D_{a}{\partial ^{2}a \over \partial x^{2}}}$

${\displaystyle {\partial i \over \partial t}=g_{a}A+g_{i}I+D_{i}{\partial ^{2}i \over \partial x^{2}}}$

Estas derivadas parciais provém das equações discretas em um domínio unidimensional para as concentrações de morfógenos quando ${\displaystyle dt}$ e ${\displaystyle dx}$ tendem à zero, de modo que a primeira parte da igualdade - uma derivada temporal de primeira ordem - passa a corresponder à uma derivada espacial de segunda ordem.^[18] Nestas equações, ${\displaystyle fa}$, ${\displaystyle fi}$, ${\displaystyle ga}$ e ${\displaystyle gi}$ são constantes e retratam interações entre o ativador e inibidor.

Generalizando para três dimensões, temos:

${\displaystyle {\partial a \over \partial t}=f(a,i)+dp\nabla ^{2}a}$

${\displaystyle {\partial i \over \partial t}=g(a,i)+dq\nabla ^{2}i}$

com ${\displaystyle \nabla ^{2}a=((\partial ^{2}/\partial x^{2})a+(\partial ^{2}/\partial y^{2})a+(\partial ^{2}/\partial z^{2})a)=\Delta a}$; e ${\displaystyle a=a(x,y,z,t)}$ e ${\displaystyle i=i(x,y,z,t)}$.

Como as equações acima lidam com a distribuição espacial de duas variáveis no tempo, a dinâmica do sistema é complexa e requer mais simplificações para um entendimento didático. Podemos aplicar uma transformada de Fourier, separando as equações de onda em seus componentes. Considerando então o número de onda ${\displaystyle k}$ (frequência angular espacial), temos que:

${\displaystyle a=a_{0}sin(kx)}$ e ${\displaystyle i=i_{0}sin(kx)}$ durante ${\displaystyle t=0}$ e que a taxa de mudança de ${\displaystyle a}$ é proporcional portanto à ${\displaystyle sin(kx)}$, conforme observado nos termos:

${\displaystyle f_{a}a_{0}sin(kx)+f_{i}i_{0}sin(kx)=(f_{a}a_{0}+f_{i}i_{0})sin(kx)}$ (Reação); e

${\displaystyle D_{a}{\partial ^{2} \over \partial x^{2}}a_{0}sin(kx)=-k^{2}D_{a}a_{0}sin(kx)}$ (Difusão).

Deste modo, podemos perceber que o número de onda ${\displaystyle k}$ de um termo não influencia o número de onda do outro, tornando-os números que evoluem de modo independente, o que facilita a análise.

Fisicamente, a evolução de número de onda é dada por:

${\displaystyle {\overrightarrow {a_{0}}}e^{\lambda t}sin(kx)}$

o que significa que um componente de onda cresce ou decai dependendo do sinal de ${\displaystyle \lambda }$. Se ${\displaystyle \lambda >0}$, temos que o componente de onda aumentará exponencialmente com o tempo. Se ${\displaystyle \lambda <0}$, o componente tende a decair de qualquer estado perturbado à um estado inicial de homogeneidade espacial.

No entanto, ${\displaystyle \lambda }$ é dependente de vários, parâmetros, como ${\displaystyle k}$, com o qual forma uma relação conhecida como Relação de Dispersão. Como a Relação de Dispersão descreve uma forma particular, na qual um arco transpassa para a região de valores positivos de ${\displaystyle \lambda }$, existe um intervalo de valores de ${\displaystyle k}$ para os quais ${\displaystyle \lambda }$ é positivo, com ${\displaystyle \lambda }$ negativo nos valores remanescentes, descrevendo uma alternância entre crescimento e caimento dos componentes de onda independentes e, portanto, das concentrações das substâncias ativadora e inibidora.^[18]