Modelo de Fitzhugh-Nagumo para o potencial de ação de neurônios

Grupo: Bernardo Boatini, Murilo Kessler Azambuja e Natália Ferrazzo

O objetivo deste trabalho é implementar e estudar a dinâmica do modelo Fitzhugh-Nagumo para potenciais de ação em células e tecidos excitáveis. O método computacional utilizado para resolver os problemas e implementar o modelo foi o FTCS (Forward Time Centered Space) e o método de Crank-Nicolson.

Potencial de Ação em Neurônios

A células vivas são sistemas eletricamente sensíveis, ou seja, podem reagir a estímulos elétricos. Isso se dá devido ao fato de que substâncias carregadas estão naturalmente vinculadas a seus processos internos de interação com o ambiente, principalmente por intermédio de canais iônicos e proteínas transmebrana como, por exemplo, a Bomba de Sódio e Potássio (Bomba Na⁺/K⁺ ATPase).

Naturalmente todas as células vivas possuem um potencial de repouso (PR) elétrico, ou seja, uma diferença de potencial elétrico, em relação ao meio (cerca de 0,1); mantida por um equilíbrio químico de concentração de íons dentro e fora da membrana plasmática.

Existem células que reagem estímulos elétricos apenas reestabelecendo o PR original por transporte passivo (sem gasto de energia) através da membrana, e estas são ditas células não-excitáveis.

Por outro lado, existem células que sob a ação do mesmo estímulo produzem um tipo de resposta bem característica: potencial de ação (PA); um pulso elétrico intenso (capaz de inverter a polarização do Potencial de Membrana) que se propaga ao longo da membrana da célula, sustentado por uma cadeia de transportes ativos (com gasto de energia) e que não decai ao longo do tempo e espaço; a esse tipo de células damos o nome de excitáveis.

Os Neurônios são as células excitáveis do tecido nervoso (que constituem o encéfalo e medula espinhal, gânglios e nervos do reino animal) e com já vimos são capazes de gerar PA. Um potencial de ação pode assumir diversos formatos, mas ao longo do axônio (Figura 1) de um neurônio eles tendem a uma curva como a da Figura 2.

Figura 1 -Representação de um potencial de ação(vermelho) ao longo de um axônio de neurônio, partindo do soma neural em direção a arvore dentrítica.

Figura 2 -Curva de um Potencial de Ação genérico no tempo, em um ponto do axônio de um neurônio.

Olhando para Figura 2 vemos alguns aspectos importantes:

• O potencial de ação necessita de um estímulo mínimo (limiar) para ser ativado, abaixo desse valor o estímulo decai como em uma célula não excitável;
• Acima desse limiar a célula segue o principio de "Tudo ou Nada", ou seja, assume o valor máximo possivel dentro de sua capacidade, independente do estímulo aplicado;
• A etapa de despolarização (crescimento) é brusca e varia mais rapidamente que a repolarização (decaimento);
• O período que contém a repolarização e hiperpolarização da membrana é chamado período refratário, e se caracteriza por não permitir que ocorra nenhum disparo até que a membrana atinja o potencial de repouso.

Modelo

Premissa do modelo

Para iniciar a modelagem do sistema, devemos antes enfatizar três condições básicas que o potencial deve obedecer para que seja um PA :

• Deve existir um limiar de voltagem para que um estímulo desencadeie o PA;
• Uma vez atingido o limiar, a voltagem deve aumentar até o máximo possível;
• Caso o estímulo não atinja o limiar, ele deve desaparecer rapidamente.

Podemos então definir $\upsilon$ como a variável normalizada que fará o papel da diferença de potencial elétrico no axônio, sendo que $\upsilon =0$ é definido como o potencial de repouso, $\upsilon =1$ é a diferença de potencial máxima suportada pela célula excitável e $\upsilon =\alpha$ é o limiar de voltagem $(0<\alpha <1)$. Desta forma podemos escrever as condições acima como

${\begin{cases}Se\ \upsilon >\alpha \Rightarrow \upsilon \rightarrow 1\\Se\ \upsilon <\alpha \Rightarrow \upsilon \rightarrow 0\\\end{cases}}\qquad (*)$

Pode-se dizer que estas condições impõem que $v=0$ e $v=1$ sejam pontos de equilíbrio estável, enquanto $v=\alpha$ seja um ponto de equilíbrio instável.

Assim, a variação temporal do potencial elétrico na célula pode ser dada por

${\frac {\partial \upsilon }{\partial t}}=f(\upsilon )\qquad (1)$

Onde $f(\upsilon )$ é alguma função que faz $v$ satisfazer as condições $(*)$.

No modelo de Nagumo, utiliza-se o polinômio de terceiro grau $f_{N}(\upsilon )=-v(v-\alpha )(v-1)$.

Substituindo $f_{N}(v)$ em $(1),$ é fácil notar que $v=0,$ $v=\alpha$ e $v=1$ são pontos de equilíbrio de $v$. Por outro lado, como podemos ver na Figura 3, se $0<v<\alpha$, temos que ${\frac {\partial \upsilon }{\partial t}}<0,$ logo o valor de $v$ diminui até atingir o ponto de equilíbrio em $v=0$. Se $\alpha <v<1,$ ${\frac {\partial \upsilon }{\partial t}}>0,$ fazendo o potencial $v$ crescer até atingir o valor máximo em $v=1$. Desta forma, temos que $v=1$ e $v=0$ são pontos de equilíbrio estável —representando os pontos de potencial máximo e mínimo, respectivamente—, enquanto que $v=\alpha$ é um ponto de equilíbrio instável, funcionando como o limiar de voltagem para desencadear o PA.

Figura 3 — Curva de $f_{N}(v)\times v$ com $\alpha =0.4$

Em outras palavras, vemos o caráter "tudo ou nada" do potencial deste modelo. Se o valor do estímulo inicial se encontrar entre $0$ e $\alpha$, o estímulo desaparece, porém se o estímulo inicial estiver entre $\alpha$ e $1$, o estímulo é amplificado até o valor máximo. Entretanto, modelando o sistema apenas com a equação $(1)$, vemos que o PA não é propagado por todo axônio, ele atua apenas localmente, onde ocorreu o estímulo. Além disso, uma vez estimulado, o neurônio nunca volta para o seu estado de repouso, permanecendo permanentemente no valor máximo de potencial.

Equação de Nagumo

A equação de Nagumo complementa o modelo acima adicionando um termo difusivo à equação $(1)$. Assim temos que a equação de Nagumo é dada por

${\frac {\partial v}{\partial t}}={\frac {\partial ^{2}v}{\partial x^{2}}}-v(v-\alpha )(v-1)\qquad (2)$

Onde $0<\alpha <{\frac {1}{2}},\ x\in \Re ,$ e $t>0$.

Desta forma o PA não é mais restrito à região onde ocorre o estímulo e se difunde ao longo de todo axônio.

Modelo de Fitzhugh-Nagumo

O modelo de Fitzhugh-Nagumo (FN) complementa a equação de Nagumo introduzindo uma nova variável, transformando a equação original em um sistema de equações. A variável $v$, chamada de "variável rápida", representa a variação do potencial na membrana da célula, enquanto que a nova variável $w$, chamada de "variável lenta" ou variável de recuperação, representa a capacidade da célula retornar ao seu PR após ser excitada por um estímulo externo. o modelo de FN é dado pelo sistema de equações

${\begin{cases}{\frac {\partial v}{\partial t}}={\frac {\partial ^{2}v}{\partial x^{2}}}-v(v-\alpha )(v-1)-w+I,\qquad \\{\frac {\partial w}{\partial t}}=\epsilon (v-\gamma w)\end{cases}}\qquad (3)$

Onde $0<\alpha <1$ é o limiar de potencial, $\epsilon$ e $\gamma$ são parâmetros a serem ajustados —relacionados à velocidade de atuação da variável de recuperação— e $I$ é a magnitude do estímulo externo .

Método de Crank-Nicolson

(explicar o método --> Natália)

Equação de Nagumo 1D

(aplicar o método na equação e testar estabilidade --> Natália)

Método FTCS

(explicar o método --> Murilo)

Equação de Recuperação 1D

(aplicar o método na equação e testar estabilidade --> Murilo))

Modelo FitzHung-Nagumo 2D

(aplicar o método na equação e testar estabilidade --> Bernardo)

O sistema de EDP's em 2 dimensões, assumindo uma difusão isotrópica, é dado por

${\frac {\partial v}{\partial t}}=D\left({\frac {\partial ^{2}v}{\partial x^{2}}}+{\frac {\partial ^{2}v}{\partial y^{2}}}\right)+f_{N}(v)-w$

${\frac {\partial w}{\partial t}}=\epsilon \left(v-\gamma w\right)$

Como a equação de recuperação já foi discretizada, e não depende da dimensão do problema, precisamos apenas aplicar o FTCS na equação difusiva de naumo em 2D, assumindo $\Delta x=\Delta y$, temos

${\frac {v_{i,j}^{n+1}-v_{i,j}^{n}}{\Delta t}}={\frac {D}{(\Delta x)^{2}}}\left[(v_{i-1,j}^{n}-2v_{i,j}^{n}+v_{i+1,j}^{n})+(v_{i,j-1}^{n}-2v_{i,j}^{n}+v_{i,j+1}^{n})\right]+f_{Ni,j}^{n}-w_{i,j}^{n}$

$v_{i,j}^{n+1}=v_{i,j}^{n}+{\frac {D\Delta t}{(\Delta x)^{2}}}\left[(v_{i-1,j}^{n}-2v_{i,j}^{n}+v_{i+1,j}^{n})+(v_{i,j-1}^{n}-2v_{i,j}^{n}+v_{i,j+1}^{n})\right]+\Delta tf_{Ni,j}^{n}-\Delta tw_{i,j}^{n}$

na qual os índices $i$ e $j$ se referem às coordenadas $x$ e $y$ respectivamente; e o índice $n$ referente ao tempo.

Vamos verificar as condições de estabilidade do problema por modos de Fourier, ou seja, substituindo $v_{i,j}^{n}=A^{n}e^{i(xq_{x}+yq_{y})}$

Resultados

Equação de Nagumo 1D

(explicar os testes e os gráficos/animações --> Natália)

Modelo FitzHung-Nagumo 1D

(explicar os testes e os gráficos/animações --> Murilo)

Modelo FitzHung-Nagumo 2D

(explicar os testes e os gráficos/animações --> Bernardo)

Discussão

(Contextualizar resultados --> Bernardo)