Exibir código-fonte para 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 FitzHungh-Nagumo, e das equações que o compõem, 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)<ref name=QUILLFELDT> https://www.ufrgs.br/mnemoforos/arquivos/potenciais2005.pdf Jorge A. Quillfeldt,"ORIGEM DOS POTENCIAIS ELÉTRICOS DAS CÉLULAS NERVOSAS" </ref>.

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<ref name=QUILLFELDT> https://www.ufrgs.br/mnemoforos/arquivos/potenciais2005.pdf Jorge A. Quillfeldt,"ORIGEM DOS POTENCIAIS ELÉTRICOS DAS CÉLULAS NERVOSAS" </ref>.

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.

[[Arquivo:PA_axonio.png|400px|left|thumb|center|'''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.]]
[[Arquivo:PA1.png|600px|rigth|thumb|center|'''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 <ref name=VIDEO> https://youtu.be/H9yxE9yrH5w </ref>:
*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 <math>\upsilon</math> como a variável normalizada que fará o papel da diferença de potencial elétrico no axônio, sendo que <math>\upsilon = 0</math> é definido como o potencial de repouso, <math>\upsilon = 1</math> é a diferença de potencial máxima suportada pela célula excitável e <math>\upsilon = \alpha</math> é o limiar de voltagem <math>(0<\alpha<1)</math>. Desta forma podemos escrever as condições acima como

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

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

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

<center><math> \frac{\partial \upsilon}{\partial t} = f(\upsilon)\qquad (1)</math></center>

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

No modelo de Nagumo, utiliza-se o polinômio de terceiro grau <math>f_N(\upsilon)= -v(v-\alpha)(v-1)</math>.

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

[[Arquivo:graph.png|400px|thumb|center|'''Figura 3''' — Curva de <math>f_N(v) \times v</math> com <math>\alpha=0.4</math>]]

Em outras palavras, vemos o caráter "tudo ou nada" do potencial deste modelo. Se o valor do estímulo inicial se encontrar entre <math>0</math> e <math>\alpha</math>, o estímulo desaparece, porém se o estímulo inicial estiver entre <math>\alpha</math> e <math>1</math>, o estímulo é amplificado até o valor máximo. Entretanto, modelando o sistema apenas com a equação <math>(1)</math>, 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 <math>(1)</math>. Assim temos que a equação de Nagumo é dada por<ref name=GARCIA> https://tamiu-ir.tdl.org/bitstream/handle/2152.4/60/GARCIA-THESIS-2015.pdf?sequence=1&isAllowed=y Gabriel Perry Natanni Garcia, "NUMERICAL SIMULATION OF THE NAGUMO EQUATION BY FINITE DIFFERENCEMETHOD" </ref><ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document </ref>

<center><math>\frac{\partial v}{\partial t} = \frac{\partial ^2 v}{\partial x^2}-v(v-\alpha)(v-1)\qquad (2)</math></center>

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

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 <math>v</math>, chamada de "variável rápida", representa a variação do potencial na membrana da célula, enquanto que a nova variável <math>w</math>, 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 <ref name=GARCIA> https://tamiu-ir.tdl.org/bitstream/handle/2152.4/60/GARCIA-THESIS-2015.pdf?sequence=1&isAllowed=y Gabriel Perry Natanni Garcia, "NUMERICAL SIMULATION OF THE NAGUMO EQUATION BY FINITE DIFFERENCEMETHOD" </ref>

<center><math>
\begin{cases}
\frac{\partial v}{\partial t} = \frac{\partial ^2 v}{\partial x^2}-v(v-\alpha)(v-1)-w, \qquad\\
\frac{\partial w}{\partial t}=\epsilon(v-\gamma r)
\end{cases}
\qquad(3)
</math></center>

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


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

<math>\frac{\partial w}{\partial t} = \epsilon\left(v - \gamma w\right) </math>


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 <math>\Delta x = \Delta y</math>, temos


<math>\frac{v_{i,j}^{n+1} - v_{i,j}^{n}}{\Delta t} = \frac{D}{(\Delta x)^2} \left[ (v_{i-1,j}^{n} - 2 v_{i,j}^{n} + v_{i+1,j}^{n}) + (v_{i,j-1}^{n} - 2 v_{i,j}^{n} + v_{i,j+1}^{n})\right] + f_{N i,j}^{n} - w_{i,j}^{n}</math>


<math>v_{i,j}^{n+1} = v_{i,j}^{n} + \frac{D \Delta t}{(\Delta x)^2} \left[ (v_{i-1,j}^{n} - 2 v_{i,j}^{n} + v_{i+1,j}^{n}) + (v_{i,j-1}^{n} - 2 v_{i,j}^{n} + v_{i,j+1}^{n})\right] + \Delta t f_{N i,j}^{n} - \Delta t w_{i,j}^{n}</math>

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

Vamos verificar as condições de estabilidade do problema por modos de Fourier, ou seja, substituindo <math>v_{i,j}^{n} = A^{n}e^{i(xq_x +yq_y)}</math>

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

==Programas==

==Referências==
<references/>
*5. [http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh, "FitzHugh-Nagumo model"]