Equação de Fitzhugh-Nagumo

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

O objetivo deste trabalho é implementar e estudar a dinâmica do modelo FitzHungh-Nagumo, e das equações que o compõem, para potenciais de ação em células e tecidos excitáveis. Os métodos computacionais 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)[1].

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[1].

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

Equação de Nagumo

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 do potencial elétrico no axônio, sendo que $\upsilon =0$ é definido como o potencial de repouso, $\upsilon =1$ é o estímulo máximo 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}}$

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 do potencial.

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, enquanto que $v=\alpha$ é um ponto de equilíbrio instável, assim como se esperava pelas condições do potencial.

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 quando utilizamos o polinômio $f_{N}(v)$. 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, neste modelo inicial 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 no valor máximo de potencial.

Equação de Nagumo difusiva

Para resolver o problema do modelo

Método de Crank-Nicolson

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

Equação de Nagumo Difusiva 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_{l,j}^{n+1}=v_{l,j}^{n}+{\frac {D\Delta t}{(\Delta x)^{2}}}\left[(v_{l-1,j}^{n}-2v_{l,j}^{n}+v_{l+1,j}^{n})+(v_{l,j-1}^{n}-2v_{l,j}^{n}+v_{l,j+1}^{n})\right]+\Delta tf_{Nl,j}^{n}-\Delta tw_{l,j}^{n}$

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

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

• substituindo $v_{l,j}^{n}=A^{n}e^{i(xq_{x}+yq_{y})}=A^{n}e^{i\beta }$
• usando $k={\frac {D\Delta t}{(\Delta x)^{2}}}$
• por simplicidade, como fizemos na seção 3.1, vamos linearizar $f_{Nl,j}^{n}\rightarrow -\alpha v_{l,j}^{n}=-\alpha A^{n}e^{i\beta }$
• por fim, desconsiderar $\Delta tw_{l,j}^{n}$

$A^{n+1}e^{i\beta }=A^{n}e^{i\beta }+kA^{n}\left[e^{i((x-\Delta x)q_{x}+yq_{y})}+e^{i((x+\Delta x)q_{x}+yq_{y})}+e^{i(xq_{x}+(y-\Delta y)q_{y})}+e^{i(xq_{x}+(y+\Delta y)q_{y})}-4e^{i\beta }\right]-\Delta t\alpha A^{n}e^{i\beta }$

$A^{n+1}e^{i\beta }=A^{n}e^{i\beta }+kA^{n}e^{i\beta }\left[e^{-i\Delta xq_{x}}+e^{i\Delta xq_{x}}+e^{-i\Delta yq_{y}}+e^{i\Delta yq_{y}}-4\right]-\Delta t\alpha A^{n}e^{i\beta }$

Dividindo os dois lados da equação por $A^{n}e^{i\beta }$ temos

${\frac {A^{n+1}}{A^{n}}}=1+k\left[e^{-i\Delta xq_{x}}+e^{i\Delta xq_{x}}+e^{-i\Delta yq_{y}}+e^{i\Delta yq_{y}}-4\right]-\Delta t\alpha$

Sabendo as seguintes identidades trigonométricas:

• $e^{-i\theta }+e^{i\theta }=2\cos(\theta )$
• $\cos(2\theta )=1-2\sin ^{2}(\theta )$

Aplicando consecutivamente a equação

${\frac {A^{n+1}}{A^{n}}}=1+2k\left[\cos(\Delta xq_{x})+\cos(\Delta yq_{y})-2\right]-\Delta t\alpha$

${\frac {A^{n+1}}{A^{n}}}=1+2k\left[(1-2\sin ^{2}(\Delta xq_{x}/2))+(1-2\sin ^{2}(\Delta yq_{y}/2))-2\right]-\Delta t\alpha$

${\frac {A^{n+1}}{A^{n}}}=1+4k\left[\sin ^{2}(\Delta xq_{x}/2)+\sin ^{2}(\Delta yq_{y}/2)\right]-\Delta t\alpha$

Sabemos que a condição de estabilidade por modos de Fourier é obtida quando $\left|{\frac {A^{n+1}}{A^{n}}}\right|<1$, portanto

$\left|1-\Delta t\alpha +4k\left[\sin ^{2}(\Delta xq_{x}/2)+\sin ^{2}(\Delta yq_{y}/2)\right]\right|<1$

Como $0<\Delta t<1$ e $0<\alpha <1$, no pior dos casos os termos de senos ao quadrado são 1, e a condição de estabilidade fica

${\frac {\alpha \Delta t}{8}}<k<{\frac {2+\alpha \Delta t}{8}}$

Resultados

Equação de Nagumo Difusiva 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)

Os resultados obtidos pela implementação da equação difusiva de Nagumo são análogos ao comportamento "tudo ou nada" de resposta em uma célula excitável, já que um estimulo acima do limiar levar ao aumento do potencial até um valor fixo. Embora esse valor não atinja o pico de ultrapassagem para todo estimulo inicial que esteja entre o limiar e a despolarização máxima da célula, essa peculiaridade é prevista pelo modelo[2].

Mas outra questão importante de se observar é a forma como diferentes larguras de estimulo inicial interferem na difusão do potencial elétrico na membrana. Um fenômeno que é verificável é que pulsos maiores difundem mais rápido sobre o axônio, já que os efeitos de inversão de campo elétrico na região do estímulo podem ser percebidos nas proximidades mesmo antes do estímulo chegar ali[1].

Quanto a implementação do Sistema FitzHung-Nagumo com potencial difusivo sobre a membrana, vemos que o potencial(nos gráficos de v por t) de ação é gerado com um perfil bem semelhante a um PA convencional de axônio. Para manter o sentido fisiológico escolhemos fazer as simulações para um $\alpha >0$ [1],que como mencionado na secção 2.3 não consegue reagir a estímulos contínuos, nem após o período refratário, ainda é muito confiável para simular estímulos únicos.

A maioria dos estudos no modelo FitzHung-Nagumo utiliza as soluções oscilatórias decorrentes $\alpha <0$, mas embora essa abordagem leve a uma gama maior de resultados para estímulo contínuo[2,4] o sentido fisiológico é prejudicado, principalmente pelo fato de não possuir um limiar[2].

Referências

• 1. Jorge A. Quillfeldt,"ORIGEM DOS POTENCIAIS ELÉTRICOS DAS CÉLULAS NERVOSAS"
• 2. Eugene M. Izhikevich and Richard FitzHugh, "FitzHugh-Nagumo model"
• 3. Gabriel Perry Natanni Garcia, "NUMERICAL SIMULATION OF THE NAGUMO EQUATION BY FINITE DIFFERENCEMETHOD"
• 4. Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia "Parameters Analysis of FitzHugh-Nagumo Model for a Reliable Simulation"