Física Computacional - Contribuições do usuário [pt-br]

Modelo de agentes de distribuição de riquezas

2021-05-28T19:01:37Z

Nataliaferrazzo:

'''Grupo: Bernardo Boatini, Murilo Kessler Azambuja e Natália Ferrazzo'''
==Introdução==
A física estatística, em particular a teoria cinética dos gases, fornece uma estrutura útil para descrever a complexidade das interações de mercado. Da mesma forma que um sistema físico composto de muitas partículas trocando energia via colisões binárias, os Modelos de Troca de Cinética consideram um conjunto de agentes econômicos interagentes que trocam de forma binária uma quantidade conservada chamada de riqueza.

Este trabalho tem como objetivo calcular a evolução temporal da distribuição de riqueza entre entre um certo numero de agentes, utilizando diferentes regras de interação e um critério para medir quantitativamente a desigualdade econômica no sistema.

==Modelo==
Vamos supor um sistema com <math>N</math> agentes, onde o agente <math>i</math> é caracterizado pela riqueza <math>w_i(t)</math> e pelo fator de aversão a risco <math>\beta_i</math> no tempo <math>t</math>. Podemos então definir uma troca de riqueza entre os agentes <math>i</math> e <math>j</math> selecionados aleatoriamente, supondo que <math>i</math> ganha uma riqueza <math>\Delta w</math> de <math>j</math> , como <ref name=BENHUR> https://arxiv.org/pdf/1904.05875.pdf CARDOSO, B. F.;GONÇALVEZ, S.; IGLESIAS, J. R.; "WEALTH DISTRIBUTION MODELS WITH REGULATIONS: DYNAMICS AND EQUILIBRIA" </ref>

<center><math>w_i(t+1) = w_i(t) + \Delta w \qquad </math></center>

<center><math>w_j(t+1) = w_j(t) - \Delta w.</math></center>

Em modelos de trocas binárias como esse, quando a regra não favorece nenhum dos lados, já foi demonstrado analíticamente que o estado final sempre leva a condensação de riqueza em um agente<ref name=MESTRADOBENHUR> https://www.lume.ufrgs.br/bitstream/handle/10183/217456/001121445.pdf?sequence=1 CARDOSO, B. F.; "A concentração de riqueza em sistemas de
trocas binárias não enviesadas
"</ref>, ou seja, máxima desigualdade. Para pensar formas de evitar a condensação, certas dinâmicas podem ser adicionadas ao problema, como, por exemplo, definir uma maior probabilidade do agente mais pobre ganhar na troca.

Para decidir quem ganha e quem perde riqueza durante a interação entre agentes, utiliza-se uma probabilidade de favorecer o agente mais pobre, evitando assim a condensação, i.e., o acúmulo de toda riqueza disponível em apenas um ou poucos agentes <ref name=BENHUR> https://arxiv.org/pdf/1904.05875.pdf CARDOSO, B. F.;GONÇALVEZ, S.; IGLESIAS, J. R.; "WEALTH DISTRIBUTION MODELS WITH REGULATIONS: DYNAMICS AND EQUILIBRIA" </ref>. Esta probabilidade é dada por <ref name=BENHUR></ref> <ref name=SCAFETTA> https://arxiv.org/pdf/cond-mat/0306579.pdf SCAFETTA, N.;WEST, B. J.; PICOZZI, S.; "A Trade-Investment Model for Distribution of
Wealth" </ref>

<center><math>p = \frac{1}{2} + f \times \frac{|w_i(t)-w_j(t)|}{w_i(t)+w_j(t)}, \qquad (1)</math></center>

onde <math>f</math> é chamado de ''fator de proteção social'', que varia de <math>0</math> —mesma probabilidade de ganho de riqueza para ambos os agentes— até <math>1/2</math> —máxima probabilidade de favorecer o agente mais pobre. Desta forma, a probabilidade do agente mais pobre ganhar a quantidade <math>\Delta w</math> em uma interação entre agentes é <math>p</math>, enquanto que a probabilidade do agente mais rico ganhar a mesma quantidade <math>\Delta w</math> é <math>1-p</math>. Além disso, vemos na equação (1) que quanto maior a desigualdade de riqueza (<math>w_i(t)-w_j(t)</math>), maior é a atuação de <math>f</math>. Isso nos mostra que o fator de proteção social é uma forma de simular a aplicação políticas sociais que favorecem a distribuição de renda na população.

Uma vez sorteado qual 2 agentes, deve-se determinar qual será a quantidade <math>\Delta w</math> a ser trocada por ambos. Existem diversas formas(regras de troca) de executar essa transferencia de riqueza (algumas delas encontram-se de forma detalhada em <ref name=MESTRADOBENHUR> https://www.lume.ufrgs.br/bitstream/handle/10183/217456/001121445.pdf?sequence=1 CARDOSO, B. F.; "A concentração de riqueza em sistemas de
trocas binárias não enviesadas
" </ref>). Porém neste trabalho focaremos apenas em dois tipos de regra: na ''regra do mínimo'' e na ''regra do perdedor'', enunciadas abaixo.

===Regra do Mínimo===

Nesta regra, temos que a quantidade de riqueza trocada entre os agentes é definida como <ref name=CAON> https://link.springer.com/content/pdf/10.1140/epjst/e2007-00072-4.pdf CAON, G.M.; GONÇALVEZ, S.; CARDOSO, B. F.; "The unfair consequences of equal opportunities: Comparing exchange models of wealth distribution" </ref>
<center><math>\Delta w = min[(1-\beta_i)w_i(t); (1-\beta_j)w_j]. \qquad (2)</math></center>

Esta regra muitas vezes também é chamada de ''regra justa'', pois a quantidade de riqueza trocada entre os agentes é a mesma, independente do ganhador, logo nenhum dos agentes é favorecido.

===Regra do Perdedor===

Neste caso, para tentar evitar condensações, temos que <math>\Delta w</math> é obtido apenas pela quantia arriscada pelo perdedor, desta forma temos <ref name = CAON></ref>

<center><math>\Delta w = (1-\beta_j)w_j(t), \qquad (3)</math></center>

lembrando que <math>j</math> é o agente perdedor. Desta forma, a quantidade de riqueza a ser trocada será sempre proporcional à fortuna de <math>j</math> (i.e., <math>w_j(t)</math>) e regulada por quanto o agente está disposto a arriscar (<math>\beta_j</math>), tornando a interação entre os agentes muito mais favorável para o perdedor.

===Coeficiente de Gini===
O Coeficiente de Gini é um índice frequentemente utilizado por economistas e organizações estatísticas para mensurar quantitativamente a desigualdade de distribuição de renda em uma determinada região. Ele é definido como <ref name = BENHUR></ref>

<center><math>G(t) = \frac{1}{2}\frac{\Sigma_{i,j}|w_i(t) - w_j(t)|}{N\Sigma_{i}w_i(t)}. \qquad(4)</math></center>

O índice de Gini varia de 0, quando todos os agentes possuem a mesma riqueza (i.e., desigualdade mínima), até 1, quando toda riqueza está concentrada em apenas um agente (i.e, desiguladade máxima). Este coeficiente é utilizado tanto para medir a desigualdade na distribuição de renda dos agentes da simulação, quanto como uma medida de dispersão, para determinar a estabilidade da distribuição de riqueza <ref name = BENHUR></ref>.

==Resultados==

Como unidade de tempo das simulações foi utilizado o MCS (Monte Carlo Step), definido como o menor número de passos necessários para que todos os agentes sejam sorteados<ref name = BENHUR></ref>

===Evolução temporal sem fator de proteção social===
Na '''Figura 1''' e na '''Figura 2''' temos a evolução temporal do coeficiente de Gini para a regra do perdedor e para a regra do mínimo, respectivamente. Estas simulações foram realizadas com um número de agentes <math>N = 10000</math> sob uma média sob 10 condições inicias diferentes(ensembles) para cada valor de <math>\beta</math>. As condições iniciais foram tais que a riqueza (<math>w_i</math>) está inicialmente distribuída no intervalo <math>[0,1]</math> seguindo uma distribuição aleatória uniforme.

Note que após um certo número de MCS, o índice de Gini tende a convergir para um valor estável, para cada valor de <math>\beta</math>. O tempo de simulação foi escolhido de forma que o equilíbrio fosse atingido. O primeiro MCS foi ignorado para melhor visualização das curvas em escala, já que o interesse é no estado estacionário.

<li style="display: inline-block;">[[Arquivo:Regradoperdedor_10ens_N10000.png|400px|left|thumb|center|'''Figura 1 -''' Evolução temporal para <math>\beta \in [0, 0.7]</math> utilizando a regra do perdedor sem fator de proteção social.]]</li>

<li style="display: inline-block;">[[Arquivo:Regradominimo2_10ens_N10000.png|400px|right|thumb|center|'''Figura 2 -''' Evolução temporal para <math>\beta \in [0, 0.7]</math> utilizando a regra do mínimo sem fator de proteção social.]]</li>

Na '''Figura 2''', vemos que o valor de estabilidade do índice de Gini sempre tende para o valor de desigualdade máxima, mesmo se mudarmos o valor de <math>\beta</math>. Desta forma, quando não inserimos um fator de proteção social no problema e assumimos que as trocas entre agentes ocorrem segundo a regra do mínimo, teremos sempre uma condensação da riqueza, levando a uma alta desigualdade econômica.

Por outro lado, na '''Figura 1''', vemos que quanto maior for o fator de aversão ao risco, menor será a desigualdade econômica e apenas quando <math>\beta = 0</math> ocorre condensação (indicando, pela equação (3), que em cada interação os agentes irão trocar toda a riqueza disponível para eles).

Este resultado (a regra do mínimo sempre gerar uma condensação enquanto que a regra do perdedor estabiliza em um valor <math>< 1</math> para <math>\beta>0</math>) ocorre pois, como o valor de beta é igual para todos os agentes, teremos que, por (3), os agentes com maior riqueza irão perder mais dinheiro na regra do perdedor do que os agentes com menor riqueza e, desta forma, a riqueza estará mais bem distribuída quando o equilíbrio for atingido. Por outro lado, como na regra do mínimo o valor trocado independe de quem perde, se forem sorteados agentes com uma diferença muito grande de riqueza teremos que, por (2), a quantidade de riqueza perdida pelo agente mais pobre será muito maior comparativamente com sua própria riqueza. Desta forma teremos que, após tempo suficiente, a regra do mínimo sempre irá tender para uma distribuição de renda desigual.

----

Nas '''Figuras 3 e 4''', temos outra evolução temporal do Índice de Gini. Nestes casos, as simulações foram realizadas com número de agentes<math>N</math> igual a <math>20000</math>(para uma obter um gini mais estável) sob uma média de 10 ensembles, assim como as anteriores. Além disso, as condições iniciais foram tais que tanto a riqueza (<math>w_i</math>) quanto o fator de aversão ao risco (<math>\beta_i</math>) foram inicialmente distribuídas no intervalo <math>[0,1]</math> seguindo uma distribuição aleatória uniforme. Desta forma podemos obter um resultado mais realista, onde cada agente tem um <math>\beta_i</math> diferente e aleatório.

<li style="display: inline-block;">[[Arquivo:regradoperdedor_beta_rand_10ens_N20000.png|400px|left|thumb|center|'''Figura 3 -''' Evolução temporal do Gini médio pela regra do perdedor e para <math>\beta</math> distribuído uniformemente no intervalo <math> [0, 1]</math> utilizando a regra do perdedor sem fator de proteção social.]]</li>

<li style="display: inline-block;">[[Arquivo:regradomin_beta_rand_10ens_N20000.png|400px|right|thumb|center|'''Figura 4 -''' Evolução temporal do Gini médio pela regra do mínimo para <math>\beta</math> distribuído uniformemente no intervalo <math> [0, 1]</math> utilizando a regra do mínimo sem fator de proteção social.]]</li>

Na '''Figura 3''', assim como na '''Figura 2''', vemos que, apesar de não ter nenhum fator de proteção social, a distribuição de riqueza estabiliza em um valor menor do que 1, exatamente por que na regra do perdedor os agentes com maior riqueza arriscam perder mais do que arriscam em ganhar.

Na '''Figura 4''', vemos a mesma tendência observada anteriormente quando a regra do mínimo foi utilizada, i.e., o equilíbrio tende para uma condensação. Isto ocorre pelo mesmo motivo anteriormente discutido: quando eventualmente ocorrer uma troca entre dois agentes com uma diferença de riqueza considerável, o agente mais pobre irá perder muito mais comparativamente a sua riqueza, devido à equação (2). Desta forma, depois de passado uma quantidade de tempo suficiente, toda riqueza vai estar concentrada em muito poucos agentes.

===Evolução temporal com fator de proteção social===

As simulações que serão apresentadas abaixo foram realizadas com um número de agentes <math>N = 10000</math> sob uma média entre 10 ensembles diferentes para cada valor de <math>f</math>. As condições iniciais em todas simulações foram tais que a riqueza (<math>w_i</math>) e a aversão ao risco (<math>\beta(i) </math>) estão distribuídas no intervalo <math>[0,1]</math> seguindo uma distribuição aleatória uniforme (<math>\beta(i) </math> permanece constante durante toda simulação).

Na '''Figura 5''' e na '''Figura 6''' temos a evolução temporal do coeficiente de Gini para a regra do perdedor(a esquerda) e para a regra do mínimo(a direita). Para cada valor de <math>f</math>, note que após um certo número de MCS o índice de Gini tende a convergir para um valor estável, se tomamos um tempo de simulação suficientemente grande (o primeiro MCS foi ignorado).

<li style="display: inline-block;">[[Arquivo:RP_f_rand.png|400px|left|thumb|center|'''Figura 5 -''' Evolução temporal da regra do perdedor para <math>\beta[i]</math> a partir de uma distribuição uniforme para diferentes fatores de proteção social.]]</li>

<li style="display: inline-block;">[[Arquivo:RM_f_rand.png|400px|right|thumb|center|'''Figura 6 -''' Evolução temporal da regra do mínimo para <math>\beta[i]</math> a partir de uma distribuição uniforme para diferentes fatores de proteção social.]]</li>

Comparando os dois gráficos é possivel ver que o Índice de Gini apresenta um período transiente antes atingir o regime de equilíbrio. O sistema atinge a o equilíbrio muito rápido quando utilizamos a regra do mínimo, enquanto a regra do perdedor demora mais, assim como é verificado em <ref name = BENHUR></ref>. Veja que para <math> f = 0 </math>, o sistema condensa para o caso da regra do mínimo, enquanto que para a regra do perdedor o sistema nunca condensa, estabilizando em <math> G \approx 0.8 </math> (para ser mais exato, em <math> G = 0.83 </math>) quando <math> f = 0 </math>, coincidindo com os resultados obtidos com a ausência do fator de proteção social nas '''Figuras 3 e 4'''.
----

Um dos resultados mais importantes extraídos das '''Figuras 5 e 6''' é ilustrado na '''Figura 7'''. Nela vemos um diagrama <math> \langle G \rangle \times f </math> , onde <math> \langle G (f) \rangle </math> é a média sobre cada ensemble do último valor calculado de <math>G</math> para ambas as regras. Nota-se que para <math>f<0,3</math> a regra do mínimo possui uma maior desigualdade na distribuição de riqueza, porém quando <math> f> 0,3 </math> o oposto ocorre, a partir do qual a regra do perdedor começa a gerar uma maior desigualdade na distribuição de riqueza do sistema, apesar da diferença entre os índices de Gini serem menos significantes do que aquelas na região de <math>f<0,3</math>.

<li style="display: inline-block;">[[Arquivo:Gxf.png|400px|thumb|center|'''Figura 6 -''' Diagrama expressando a relação do Índice de Gini médio para as duas regras de troca diferentes, em função do fator de proteção social.]]</li>

Esse resultado é amplamente discutido na literatura, e indica que mesmo uma troca justa como a regra do mínimo, que tende sempre a condensação, pode gerar distribuições mais igualitárias que uma troca que naturalmente favorece os agentes mais pobres (regra do perdedor). Para <math> f<0,3 </math>, a regra do mínimo sempre gera uma distribuição de riqueza mais desigual do que aquelas geradas utilizando a regra do perdedor, porém em <math> f \approx 0,3</math> as duas linhas se cruzam, ou seja, as duas regras geram a mesma distribuição de riqueza. A partir de então, para <math> f>0,3 </math>, a regra do mínimo começa a gerar distribuições mais igualitárias do que a regra do perdedor. Desta forma, temos que para valores altos de <math> f </math>, a regra do mínimo auxilia a ter uma distribuição mais "humanitária", uma vez que já exista políticas sociais suficientes que promovam uma distribuição de renda igualitária. Enquanto isso, a regra do perdedor auxilia a resgatar mais agentes da região de baixa riqueza quando o valor de <math> f </math> é muito pequeno <ref name = CAON> </ref>.

==Discussão==

Os modelos de distribuição de riqueza por trocas binarias, se mostram muito eficientes para descrever o fenômeno da desigualdade de renda na vida real. As regras do perdedor e do mínimo são propostas muito ricas, mas existem outros modelos conhecidos e passíveis de serem testados. Os resultados que apresentamos aqui são apenas algumas potencialidades que o modelo pode apresentar e condizem com os resultados da literatura consultada. Para aprimoramentos futuros poderia ser interessantes:

* Aumentar o Tempo de simulação para observar as variações no estado transiente do Gini médio, principalmente nos modelos com o fator social;
* Estudar as influências da implementação de diferentes regras de troca entre agentes, que não sejam as regras do mínimo e do perdedor <ref name = MESTRADOBENHUR> </ref>;
* Comparar as distribuições de riquezas geradas pela simulação com distribuições reais de renda per capita, por exemplo<ref name=IGLESIAS> https://arxiv.org/pdf/2005.06106.pdf IGLESIAS, J. R.; CARDOSO, B. F.;GONÇALVEZ, S.; "Inequality, a scourge of the XXI century" </ref>;
* Estudar como perturbações de <math>f</math> no tempo podem interferir no Gini estacionário do sistema <ref name = BENHUR> </ref>;
* Implementar medidas de mobilidade financeira (chamada de liquidez) e correlação de riqueza, que só podem ser calculadas entre MCS <ref name = BENHUR> </ref> <ref name = CAON> </ref> <ref name = MESTRADOBENHUR> </ref>;
* Testar outras formas mais específicas de auxiliar os agentes mais pobres, que não seja apenas pelo fator de proteção social, como por exemplo taxação uniforme e taxação proporcional à riqueza<ref name = IGLESIAS> </ref>.

==Programas==

<source lang="python">

import numpy as np
import matplotlib.pyplot as plt

def Gini(W):
N = len(W)
I = 2 * np.arange(1, N + 1) - N - 1
return np.sum(I * np.sort(W)) / (N * np.sum(W))

</source>

===Regra do Perdedor sem Fator de Proteção Social===

<source lang="python">

'''parametros'''
N = 10000
n = 1000
GF = 0
e = 10
Gm = np.zeros((e, n))
a_array = np.arange(0, 1, 0.1)
for a in a_array:
for v in range(e):
'''microestado'''

W = np.random.rand(N)

'''macroestado'''
Go = Gini(W)
G = np.zeros(n)

for t in range(n * N):

i = np.random.randint(N)
j = np.random.randint(N)

delta_w = (1-a) * W[j]

W[i] = W[i] + delta_w
W[j] = W[j] - delta_w

if (t % N == 0):
G[t // N] = Gini(W)

T = np.arange(n)
GF = GF + Gini(W)
Gm[v] = G

GM = (GF/e)*np.ones(n)
Gm = np.sum(Gm, axis = 0)/e

plt.plot(T,Gm, label = "$\beta$ = " + str("%.2f" % a) + "")
print("para beta = " + str(a) + " GM = " + str("%.2f" % GM[n-1]) + "")
plt.legend(prop={'size': 8})
plt.xlabel("MCS")
plt.ylabel("$\langle G \rangle$")
plt.ylim(0, 1)
plt.xlim(left = 1)
Gm = np.zeros((e, n))
GF = 0
GM = 0

plt.show()

</source>

===Regra do Mínimo sem Fator de Proteção Social===

<source lang="python">

'''parametros'''
N = 10000
n = 1000
GF = 0
e = 10
Gm = np.zeros((e, n))
a_array = np.arange(0, 1, 0.1)
for a in a_array:
for v in range(e):
'''microestado'''

W = np.random.rand(N)

'''macroestado'''
Go = Gini(W)
G = np.zeros(n)

for t in range(n * N):

i = np.random.randint(N)
j = np.random.randint(N)

delta_w = (1-a) * min(W[i], W[j])

W[i] = W[i] + delta_w
W[j] = W[j] - delta_w

if (t % N == 0):
G[t // N] = Gini(W)

T = np.arange(n)
GF = GF + Gini(W)
Gm[v] = G

GM = (GF/e)*np.ones(n)
Gm = np.sum(Gm, axis = 0)/e

plt.plot(T,Gm, label = "$\beta$ = " + str("%.2f" % a) + "")
print("para beta = " + str(a) + " GM = " + str("%.2f" % GM[n-1]) + "")
plt.legend(prop={'size': 8})
plt.xlabel("MCS")
plt.ylabel("$\langle G \rangle$")
plt.ylim(0, 1)
plt.xlim(left = 1)
Gm = np.zeros((e, n))
GF = 0
GM = 0

plt.show()

</source>

===Regra do Perdedor com Fator de Proteção Social===
<source lang="python">

'''parametros'''
N = 10000 # Numero de Agentes
n = 10000 # Numero de MCS
GF = 0
e = 10 # Numero de amostras
a = np.random.rand(N) # Aversão ao risco
Gm = np.zeros((e, n))
f_array = np.arange(0, 0.6, 0.1) # Fatores de proteção social

for f in f_array:
for v in range(e):
'''microestado'''

W = np.random.rand(N) # Riqueza inicial dos agentes

'''macroestado'''
Go = Gini(W)
G = np.zeros(n)

for t in range(n * N):

i = np.random.randint(N)
j = np.random.randint(N)

# Regra do Perdedor com Fator de Proteção Social
p = 0.5+f*((abs(W[i]-W[j]))/(W[i]+W[j]))
aux = np.random.rand()

if aux < p:
if W[j]>W[i]:
delta_w = (1 - a[j]) * W[j]
W[i] = W[i] + delta_w
W[j] = W[j] - delta_w
else:
delta_w = (1 - a[i]) * W[i]
W[i] = W[i] - delta_w
W[j] = W[j] + delta_w
else:
if W[j]<W[i]:
delta_w = (1 - a[j]) * W[j]
W[i] = W[i] + delta_w
W[j] = W[j] - delta_w
else:
delta_w = (1 - a[i]) * W[i]
W[i] = W[i] - delta_w
W[j] = W[j] + delta_w

if (t % N == 0):
G[t // N] = Gini(W)

T = np.arange(n)
GF = GF + Gini(W)
Gm[v] = G

GM = (GF/e)*np.ones(n)
Gm = np.sum(Gm, axis = 0)/e

# Gráficos
plt.plot(T,Gm, label = "f = " + str("%.2f" % f) + "")
print("para f = " + str(f) + " --> GM = " + str("%.2f" % GM[n-1]) + "") # imprime os valores para o gráfico de G x f
plt.legend(prop={'size': 8})
plt.xlabel("MCS")
plt.ylabel("$\\langle G \\rangle$")
plt.xlim(left = 1)
plt.title("Regra do Perdedor $\\beta$ $\\in$ [0,1]")
Gm = np.zeros((e, n))
GF = 0
GM = 0

plt.show()

</source>

===Regra do Mínimo com Fator de Proteção Social===

<source lang="python">

'''parametros'''
N = 10000 # Numero de Agentes
n = 10000 # Numero de MCS
GF = 1
e = 10 # Numero de amostras
a = np.random.rand(N) # Fatoes de Aversão ao Risco
Gm = np.zeros((e, n))
f_array = np.arange(0, 0.6, 0.1) # Fatores de Proteção Social

for f in f_array:
for v in range(e):
'''microestado'''

W = np.random.rand(N) # Riqueza inicial dos agentes

'''macroestado'''
Go = Gini(W)
G = np.zeros(n)

for t in range(n * N):

i = np.random.randint(N)
j = np.random.randint(N)

# Regra do MÍnimo com Fator de Proteção Social
p = 0.5+f*((abs(W[i]-W[j]))/(W[i]+W[j]))
aux = np.random.rand()
delta_w = min((1 - a[i]) * W[i], (1 - a[j]) * W[j])

if aux < p:
if W[j]>W[i]:
W[i] = W[i] + delta_w
W[j] = W[j] - delta_w
else:
W[i] = W[i] - delta_w
W[j] = W[j] + delta_w
else:
if W[j]<W[i]:
W[i] = W[i] + delta_w
W[j] = W[j] - delta_w
else:
W[i] = W[i] - delta_w
W[j] = W[j] + delta_w

if (t % N == 0):
G[t // N] = Gini(W)

T = np.arange(n)
GF = GF + Gini(W)
Gm[v] = G

GM = (GF/e)*np.ones(n)
Gm = np.sum(Gm, axis = 0)/e

# Gráficos
plt.plot(T,Gm, label = "f = " + str("%.2f" % f) + "")
print("para f = " + str(f) + " --> GM = " + str("%.2f" % GM[n-1]) + "") # imprime os valores para o gráfico de G x f
plt.legend(prop={'size': 8})
plt.xlabel("MCS")
plt.ylabel("$\\langle G \\rangle$")
plt.xlim(left = 1)
plt.title("Regra do Mínimo $\\beta$$\\in$[0,1]")
Gm = np.zeros((e, n))
GF = 0
GM = 0

plt.show()

</source>

==Referências==
<references/>

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

2021-04-05T14:24:32Z

Nataliaferrazzo:

'''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)<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| A simple spiking neuron model: sodium channels alone </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 Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </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>
<ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </ref>
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>

<center><math>
\begin{cases}
\frac{\partial v}{\partial t} = D\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)
</math></center>

Onde <math>0<\alpha<1</math> é o limiar de potencial, <math>D</math> é uma constante de difusão, <math>\epsilon</math> e <math>\gamma</math> são parâmetros positivos a serem ajustados —relacionados à velocidade de atuação da variável de recuperação— e <math>I</math> é a magnitude do estímulo externo <ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </ref>
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>

== Método Crank-Nicolson ==

Uma equação diferencial parcial da forma
<center><math>Au_{xx} + Bu_{xy} + Cu_{yy} = f(x, y, u_x, u_y) \qquad (3.1)</math></center>

onde A, B e C são constantes, é chamada de equação quasilinear.
Um método implícito no tempo e estável numericamente de solucionar este tipo de equação foi proposto na metade do século XX por John Crank e Phyllis Nicolson.
O método de Crank-Nicolson é de segunda ordem no tempo e no espaço e é baseado em diferenças centradas no espaço e regra trapezoidal no tempo.
Como provaremos mais adiante, este método é incondicionalmente estável.

A equação diferencial finita de Nagumo é dada pela segunda derivada espacial obtida a partir de uma combinação da derivada temporal no passo <math> n</math> e <math> n+1</math>:

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} = (\theta) \frac{u_{i+1}^n -2u_i^{n} + u_{i+1}^n}{\Delta x^2} + (1-\theta) \frac{u_{i+1}^{n+1} -2u_i^{n+1} + u_{i-1}^{n+1}}{\Delta x^2} \qquad(3.2)</math></center>

A partir da equação (3.2), podemos obter os diferentes métodos de diferenciação da equação de Nagumo ajustando o valor de <math>\theta</math>.
*<math>\theta = 1</math> é o método explícito.
*<math>\theta = 0</math> é o método implícito.
*<math>\theta = \frac{1}{2}</math> é o método Crank-Nicolson.

Substituindo o valor de <math>\theta</math> na equação (3.2), temos:

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} = \frac{u_{i+1} ^n -2u_i ^{n} + u_{i+1} ^n}{\Delta x^2} + u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.3)</math></center>

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} = \frac{u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1}}{\Delta x^2} + u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha) \qquad(3.4)</math></center>

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} =\frac{u_{i+1} ^n -2u_i ^{n} + u_{i+1} ^n}{2\Delta x^2} + \frac{u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1}}{2\Delta x^2} + u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.5)</math></center>

Para nosso problema, usaremos o método pseudo Crank-Nicolson, devido ao termo não-linear presente em (3.5).

=== Análise de Estabilidade de Von Neumann ===
Baseada nas deduções de G. Garcia <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>.

Ao analisar a estabilidade do método, precisamos verificar se o erro gerado no procedimento não aumenta quando testamos os extremos da solução.
Para realizar esta verificação, vamos utilizar o método de Von Neumann.
A primeira coisa a se fazer é linearizar o termo não linear da equação (3.5):

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} =\frac{u_{i+1} ^n -2u_i ^{n} + u_{i+1} ^n}{2\Delta x^2} + \frac{u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1}}{2\Delta x^2} -\alpha u_{i} ^{n}\qquad(3.6)</math></center>

Se substituírmos <math>u_{i} ^{n} = \Psi (t)e^{j\gamma \Delta x}</math></div> no tempo <math>t</math>, onde <math>j = \sqrt{(-1)}</math> e <math>\gamma > 0</math> o método se mostra estável pelo seguinte:

<center><math>\bigg| \frac{\Psi (t + \Delta t)}{\Psi (t)} \bigg| \le 1 </math></center>

Substituindo <math>u_{i} ^{n} = \Psi (t)e^{j\gamma x} </math> na equação (3.6), temos:

<center><math>\begin{matrix} \frac{\Psi (t + \Delta t) e^{j\gamma x}}{\Delta t}& =&\frac{\Psi (t) e^{j\gamma (x + \Delta x)} -2 \Psi (t) e^{j\gamma x} + \Psi (t) e^{j\gamma (x - \Delta x)}}{2\Delta x^2} \\ \ &+ &\frac{\Psi (t + \Delta t) e^{j\gamma (x + \Delta x)} -2\Psi (t + \Delta t)e^{j\gamma x} + \Psi (t + \Delta t) e^{j\gamma (x + \Delta x)}}{2\Delta x^2} -\alpha \Psi (t) e^{j\gamma x}\end{matrix}\qquad(3.7)</math></center>

Substituindo <math>r = {\Delta t}/{\Delta x^2}</math> em (3.7), temos:

<center><math>\begin{matrix} \Psi (t + \Delta t) e^{j\gamma x} &=& r \bigg(\frac{\Psi (t) e^{j\gamma (x + \Delta x)} -2\Psi (t) e^{j\gamma x} + \Psi (t) e^{j\gamma (x - \Delta x)}}{2}\bigg) \\ \ &
+& r\bigg( \frac{\Psi (t + \Delta t) e^{j\gamma (x + \Delta x)} -2\Psi (t + \Delta t)e^{j\gamma x} + \Psi (t + \Delta t) e^{j\gamma (x + \Delta x)}}{2}\bigg) -\Delta t \alpha \Psi (t) e^{j\gamma x}\end{matrix}\qquad(3.8)</math></center>

Podemos notar que há <math>e^{j\gamma x} </math> em todos os termos. Simplificando-os e movendo todos os termos com <math>\Psi (t + \Delta t) </math> para o lado esquerdo e todos com <math>\Psi (t) </math> para a direita, a equação (3.8) toma a forma:

<center><math>\Psi (t + \Delta t) - r \bigg(\frac{\Psi (t + \Delta t)e^{j\gamma \Delta x} -2\Psi (t + \Delta t) + \Psi (t + \Delta t)e^{-j\gamma \Delta x}}{2}\bigg) = r \bigg( \frac{\Psi (t) e^{j\gamma \Delta x} -2\Psi (t) + \Psi (t) e^{-j\gamma \Delta x}}{2}\bigg) + \Psi (t) -\Delta t \alpha \Psi (t)\qquad(3.9)</math></center>

Fatorando os termos, ficamos com:

<center><math>\Psi (t + \Delta t) \bigg( 1 - r \bigg( \frac{e^{j\gamma \Delta x} -2 + e^{-j\gamma \Delta x}}{2}\bigg) \bigg) = \Psi (t) \bigg( r \bigg( \frac{e^{j\gamma \Delta x} - 2 + e^{-j\gamma \Delta x}}{2} \bigg) + 1 -\Delta t \alpha \bigg) \qquad(3.10)</math></center>

Utilizando a identidade trigonométrica <math> cos(x) = \frac{e^{jx} + e^{-jx}}{2}</math>, ficamos com:

<center><math>\Psi (t + \Delta t) (1 - r (cos(\gamma \Delta x) -1)) = \Psi (t) (r (cos(j\gamma \Delta x) -1 ) +1 -\Delta t \alpha)\qquad(3.11)</math></center>

Aplicando <math>1 - cos(2u) = 2 sin^2 (u)</math>, temos:

<center><math>\Psi (t + \Delta t) \bigg(1 + 2r sin^2 \bigg( \frac{\gamma \Delta x}{2}\bigg) \bigg) = \Psi (t) \bigg(1 - 2r sin^2 \bigg( \frac{\gamma \Delta x}{2}\bigg) -\Delta t \alpha\bigg)\qquad(3.12)</math></center>

Pode-se demonstrar que, a partir da equação (3.12), tem-se:

<center><math>\frac{\Psi (t + \Delta t) e^{j\gamma x}}{\Delta t} =\frac{1 - 2r sin^2( \frac{\gamma \Delta x}{2}) -\Delta t \alpha}{1 + 2r sin^2( \frac{\gamma \Delta x}{2})}\qquad(3.13)</math></center>

Como:

<center><math>\bigg| \frac{1 - 2r sin^2( \frac{\gamma \Delta x}{2}) -\Delta t \alpha}{1 + 2r sin^2( \frac{\gamma \Delta x}{2})}\bigg| \le 1 \qquad(3.13)</math></center>

prova-se que o método Crank-Nicolson é estável para <math>\Delta t </math> suficientemente pequenos.

=== Implementação para a equação de FitzHung-Nagumo ===

Para a implementação na equação de FitzHung-Nagumo, as seguintes derivadas parciais são utilizadas:

<center><math>u(x_{i} , t_{n}) = u_{i}^n \qquad(3.14)</math></center>

<center><math>\frac{\partial u}{\partial t} (x_{i} , t_{n})= \frac{u_{i}^{n+1}- u_{i}^{n}}{\Delta t}\qquad(3.15)</math></center>

<center><math>\frac{\partial^2 u}{\partial^2 x} (x_{i} , t_{n})= \frac{1}{2}\bigg(\frac{u_{i+1}^{n+1}- 2u_{i}^{n+1}+u_{i-1}^{n+1}}{\Delta x^2} + \frac{u_{i+1}^{n}- 2u_{i}^{n}+u_{i-1}^{n}}{\Delta x^2}\bigg)\qquad(3.16)</math></center>

Neste ponto, utilizamos uma aproximação para a derivada temporal e a derivada espacial de segunda ordem e substituímos as equações (3.14), (3.15) e (3.16) na equação (2).

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} =\frac{1}{2\Delta x^2}(u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1})+ \frac{1}{2\Delta x^2}(u_{i+1} ^{n} -2u_i ^{n} + u_{i-1} ^{n})
u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.17)</math></center>

Multiplicando por <math>\Delta t</math>, temos:

<center><math>u_i ^{n+1} - u_i ^{n} =\frac{\Delta t}{2\Delta x^2}(u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1})+ \frac{\Delta t}{2\Delta x^2}(u_{i+1} ^{n} -2u_i ^{n} + u_{i-1} ^{n})
\Delta t u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.18)</math></center>

Fazendo a substituição:

<center><math> r = \frac{2\Delta t}{\Delta x^2}</math></center>

e rearranjando, temos:

<center><math>-ru_{i+1}^{n+1} + (1+2r)u_{i}^{n+1} - ru_{i-1}^{n+1} = r(u_{i+1}^{n} + u_{i-1}^{n}) + (1-2r)u_{i}^{n} + \Delta tu_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.19)</math></center>

Podemos, então representar a equação (3.19) em forma matricial:

<center><math>\mathbf{AU} = \mathbf{F}</math></center>

Onde o lado esquerdo da equação é a matriz tridimensional a seguir
<center><math>
A =
\begin{bmatrix}
1+2r & -r & 0 & \cdots & 0 \\
-r & 1+2r & -r & \cdots & 0 \\
0 & -r & 1+2r & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & -r \\
0 & 0 & 0 & -r & 1+2r \\
\end{bmatrix}
</math></center>

e o lado direito:
<center><math>
F =
\begin{bmatrix}
0 \\
\vdots \\
r(u_{n+1 ^n} + u_{i-1 ^n}) + (1-2r)u_i ^n + \Delta t u_i ^n(1-u_i ^n)(u_i ^n - \alpha) \\
\vdots \\
1\\
\end{bmatrix}
</math></center>

Com isso, podemos encontrar as soluções da equação resolvendo:

<center><math>\mathbf{U} = \mathbf{F/A}</math></center>

===FTCS Explícito===
No método explícito, a discretização é

<math>
\begin{cases}
\frac{\partial f}{\partial t} \rightarrow \frac{f(x, t+\Delta t) - f(x,t)}{\Delta t}\\
\frac{\partial^2 f}{\partial x^2} \rightarrow \frac{f(x, t+\Delta t) - 2f(x,t) + f(x-\Delta t}{(\Delta x)^2}
\end{cases}
</math>

Substituindo na EDU e utilizando a notação de índices (<math>n</math> é índice temporal e <math>j</math> é o índice temporal), temos

<math>
\frac{f_j^{n+1} - f_j^n}{\Delta t} = D\frac{f_{j+1}^n - 2f_j^n + f_{j-1}^n}{(\Delta x)^2}
</math>

Agrupando <math>\Delta t,\ \Delta x\ e\ D</math> em uma constante <math>k</math> e deixando o passo temporal futuro em função do passo anterior explicitamente, temos

<math>
f_j^{n+1} = f_j^n+k(f_{j+1}^n - 2f_j^n + f_{j-1}^n).
</math>

Utilizando esta equação é possível calcular o valor de <math>f</math> iterativamente para todo intervalo de tempo.

===FTCS Implícito===
No método implícito, também chamado de BTCS (para trás no tempo centrado no espaço, do inglês: Backwards Time Centered Space), pois o passo temporal é dado para trás, assim a derivada temporal é

<math>
\begin{cases}
\frac{\partial f}{\partial t} \rightarrow \frac{f_j^{n} - f_j^{n-1}}{\Delta t}\\
\end{cases}
</math>

Resultando na seguinte EDU

<math>
f_j^{n} = f_j^{n-1}+k(f_{j+1}^n - 2f_j^n + f_{j-1}^n).
</math>

Reescrevendo os índices <math>n</math> temos

<math>
f_j^{n+1} = f_j^{n}+k(f_{j+1}^{n+1} - 2f_j^{n+1} + f_{j-1}^{n+1})
</math>

Onde <math>f_j^{n+1}</math> está escrita implicitamente, e pode ser calculada utilizando o Algoritmo de Thomas [https://moodle.ufrgs.br/pluginfile.php/3653465/mod_resource/content/1/thomas.pdf].

===Equação de Recuperação 1D===

No modelo de Fitzhugh-Nagumo, utilizou-se o mesmo algoritmo de Crank-Nicolson utilizado na Equação de Nagumo, na seção 3.1. Entretanto, ainda é necessário resolver a equação diferencial de <math>w</math>, dada por

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

Para isso utilizou-se o método FTCS explícito. Discretizando a derivada, temos

<math>
\frac{\partial w}{\partial t} \rightarrow \frac{w(x,t + \Delta t)-w(x,t)}{\Delta t}
</math>

Utilizando a notação de índices e aplicando na equação da recuperação, temos

<math>
\frac{w_j^{n+1}-w_j^n}{\Delta t} = \epsilon(v_j^n - \gamma w_j^n)
</math>

<math>
w_j^{n+1}-w_j^n = \Delta t \epsilon(v_j^n - \gamma w_j^n).
</math>

Assim, uma vez conhecido <math>v_j^n</math>, é possível calcular <math>w</math> para qualquer intervalo de tempo.

Vamos verificar as condições de estabilidade do método utilizando o os modos de Fourier:
*substituindo <math>w_j^n=A^ne^{i(jq\Delta x)},</math> onde <math>i</math> é a unidade imaginária;
*desconsiderando <math>v_j^n</math>,
temos

<math>
A^{n+1}e^{i(jq \Delta x)} = -\Delta t \epsilon \gamma e^{i(lq \delta x)}
</math>

Dividindo ambos os lados por <math>A^ne^{i(jq\Delta x)},\ </math> temos

<math>
\frac{A^{n+1}}{A^n} = -\Delta t \epsilon \gamma
</math>

Para que o método seja estável, temos que <math>\Bigg|\frac{A^{n+1}}{A^n}\Bigg| < 1,\ </math> logo

<math>
\Bigg|\frac{A^{n+1}}{A^n}\Bigg| =\Bigg|-\Delta t \epsilon \gamma \Bigg|< 1.
</math>

Assim, temos que a condição de estabilidade sobre o método é tal que

<math>
\frac{-1}{\epsilon \gamma}<\Delta t < \frac{1}{\epsilon \gamma}.
</math>

===Modelo FitzHung-Nagumo 2D===
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 de Nagumo em 2D. Assumindo <math>\Delta x = \Delta y</math>, temos

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

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

<math>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} - 4 v_{l,j}^{n}\right] + \Delta t f_{N l,j}^{n} - \Delta t w_{l,j}^{n}</math>

onde os índices <math>l</math> e <math>j</math> se referem às coordenadas <math>x</math> e <math>y</math>, respectivamente, e o índice <math>n</math> refere-se ao tempo.

Vamos verificar as condições de estabilidade do método utilizando os modos de Fourier:
*substituindo <math>v_{l,j}^{n} = A^{n}e^{i(xq_x +yq_y)} = A^{n}e^{i\beta}</math>
*usando <math>r = \frac{D \Delta t}{(\Delta x)^2}</math>
*por simplicidade, como fizemos na seção 3.1, vamos linearizar <math>f_{N l,j}^{n} \rightarrow -\alpha v_{l,j}^{n} =-\alpha A^{n}e^{i\beta}</math>
*por fim, desconsiderando <math>\Delta t w_{l,j}^{n}</math>, temos

<math>A^{n+1}e^{i\beta} = A^{n}e^{i\beta} + rA^{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)} - 4 e^{i\beta}\right] -\Delta t \alpha A^{n}e^{i\beta}</math>

<math>A^{n+1}e^{i\beta} = A^{n}e^{i\beta} + rA^{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}</math>

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

<math>\frac{A^{n+1}}{A^{n}} = 1 + r\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</math>

Sabendo as seguintes identidades trigonométricas:

*<math>e^{-i\theta} + e^{i\theta} = 2\cos(\theta)</math>
*<math>\cos(2\theta) = 1 - 2\sin^2(\theta)</math>

e aplicando à equação, obtemos

<math>\frac{A^{n+1}}{A^{n}} = 1 + 2r\left[ \cos(\Delta xq_x) + \cos(\Delta yq_y) - 2\right] -\Delta t \alpha</math>

<math>\frac{A^{n+1}}{A^{n}} = 1 + 2r\left[ (1 - 2\sin^2(\Delta xq_x/2)) + (1 - 2\sin^2(\Delta yq_y/2)) - 2\right] -\Delta t \alpha</math>

<math>\frac{A^{n+1}}{A^{n}} = 1 + 4r\left[ \sin^2(\Delta xq_x/2) + \sin^2(\Delta yq_y/2)\right] -\Delta t \alpha.</math>

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

<math>\left|1 -\Delta t \alpha + 4r\left[ \sin^2(\Delta xq_x/2) + \sin^2(\Delta yq_y/2)\right]\right| < 1</math>.

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

<math>\frac{\alpha \Delta t}{8} < r < \frac{2 + \alpha \Delta t}{8}</math>.

==Resultados==
===Equação de Nagumo 1D===

Para as simulações numéricas foram utilizados os seguintes parâmetros:

* <math>\alpha = 0,1</math>
* <math>\Delta t = 0,01\qquad para 0< t < 20000</math>
* <math>\Delta x = 0,01\qquad para 0< x < 200</math>
* <math>\Delta t = 0,01\qquad para 0<\Delta t < 20000</math>

Com isso, dividimos as simulações em dois casos.

'''Caso I:''' As simulações foram geradas a partir da definição de uma largura de estímulo estreita <math>l= 20</math>.
Com isso, variamos a intensidade do estímulo. Nas seguintes figuras, podemos ver as simulações para uma intensidade menor do que o limiar, igual ao limiar e maior que o limiar, respectivamente, <math>I < \alpha</math>, <math>I = \alpha</math> e <math>I > \alpha</math>.

[[Arquivo:N_fino_abaixolimiar.gif|800px|thumb|center|'''Figura 4''' — estímulo <math>\ell = 20</math> <math>I = 0.05</math>]]

Como pode ser visto na Figura 4, quando <math>I=0,05</math>, ou seja, menor que o limiar, o sinal rapidamente se esvai.

Para uma intensidade igual ao limiar, o sinal também se extingue, porém com um intervalo de tempo maior do que anteriormente.
[[Arquivo:N_fino_iguallimiar.gif|800px|thumb|center|'''Figura 6''' — estímulo <math>\ell = 20</math> <math>I = 0.1</math>]]

E, finalmente, quando a intensidade é maior que o limiar, não importa quão maior seja, o sinal é passado ao máximo da função.
[[Arquivo:N_fino_acimalimiar.gif|800px|thumb|center|'''Figura 5''' — estímulo <math>\ell = 20</math> <math>I = 0.2</math>]]

'''Caso II:''' As seguintes simulações foram geradas a partir da definição de uma largura de estímulo larga <math>l= 100</math>.
Com isso, novamente a intensidade do estímulo é variada.

[[Arquivo:N_largo_abaixolimiar.gif|800px|thumb|center|'''Figura 7''' — estímulo <math>\ell = 100</math> <math>I = 0.05</math>]]
Para uma intensidade menor que o limiar, novamente o sinal se esvai.

[[Arquivo:N_largo_iguallimiar.gif|800px|thumb|center|'''Figura 9''' — estímulo <math>\ell = 100</math> <math>I = 0.1</math>]]
O mesmo acontece quando a intensidade é igual ao limiar.

Já quando a intensidade é maior que o limiar, o comportamento “tudo ou nada” novamente se mostra presente.
[[Arquivo:N_largo_acimalimiar.gif|800px|thumb|center|'''Figura 8''' — estímulo <math>\ell = 100</math> <math>I = 0.2</math>]]

Podemos notar que quando a largura é maior, o sinal se esvai mais lentamente se comparado a uma largura de sinal mais estreita.

===Modelo FitzHung-Nagumo 1D===
Para a implementação da simulação do modelo FN em 1D utilizou-se <math>\Delta t = 0.01</math> (com <math>0<t<30000</math>), <math>\Delta x = 1</math> (com <math>0<x<120</math>) e os demais parâmetros foram idênticos àqueles encontrados na ref.[1] <ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </ref>. Eles são

<math>
\begin{cases}
\alpha = 0.1\\
\epsilon = 0.005\\
\gamma = 0.0001\\
D = 0.1\\
\end{cases}
</math>

Foram aplicadas condições de contorno nulas nas bordas do intervalo de <math>x</math>.

[[Arquivo:FN_vXt.png|800px|thumb|center|'''Figura 10''' — Evolução temporal do potencial elétrico em um ponto da membrana do axônio quando excitado por um estímulo acima do limiar de potencial.]]

Na figura 10 vemos a evolução temporal do potencial elétrico em um ponto da membrana do axônio (<math>x=15</math>) quando excitado acima do limiar de potencial. Note que a curva do PA é muito semelhante àquela da na Figura 2, apresentando vários dos aspectos característicos, como a assimetria entre os períodos de polarização (ascensão rápida) e repolarização (decaimento mais lento) e um período refratário, onde o potencial lentamente se recupera em direção do PR. Podemos assim entender a escolha dos valores dos parâmetros <math>\epsilon</math> e <math>\gamma</math>, que controlam a velocidade de crescimento da variável de recuperação (assim como se pode ver na equação (3)). Experimentalmente, sabemos que o período refratário é muito mais longo do que o período de despolarização
<ref name=QUILLFELDT> https://www.ufrgs.br/mnemoforos/arquivos/potenciais2005.pdf Jorge A. Quillfeldt,"ORIGEM DOS POTENCIAIS ELÉTRICOS DAS CÉLULAS NERVOSAS" </ref>
, logo precisamos que a intensidade da variável <math>w</math> aumente muito lentamente. Para isso, escolhe-se valores pequenos de <math>\epsilon</math> e <math>\gamma</math>, pois caso contrário o potencial da membrana celular não teria tempo de atingir o pico do PA.

A velocidade de crescimento da variável de recuperação também explica a ausência do "tudo ou nada" no PA da figura 10. Na figura 11, vê-se a evolução temporal de um PA com <math>\epsilon = 0.001</math> (5 vezes menor do que aquele da figura 10).

[[Arquivo:FN_vXt_Epequeno.png|800px|thumb|center|'''Figura 11''' — Evolução temporal do potencial elétrico em um ponto da membrana do axônio para <math>\epsilon = 0.001</math>.]]

Na figura 11, como a variável de recuperação cresce muito mais lentamente, além de demorar mais tempo para a membrana celular sair do período refratário, vemos que o pico do PA de ação é muito mais próximo do valor máximo. Isso ocorre porque a variável de recuperação cresce mais lentamente, possibilitando que a variável rápida atinja um valor maior antes de começar a diminuir devido à ação de <math>w</math>.

[[Arquivo:FN_abaixolim_apos.gif|800px|thumb|center|'''Figura 12''' — Simulação de um estímulo excitando um axônio, segundo o modelo de Fitzhugh-Nagumo. Essa simulação foi feita para <math>\alpha = 0.1</math> e um estímulo inicial de largura <math>\ell = 10</math> e intensidade <math>I=0.05</math> (abaixo do limiar). Note que o estímulo não desencadeia um PA e logo desaparece.]]

[[Arquivo:FN_acimalim_apos.gif|800px|thumb|center|'''Figura 13''' — Simulação de um estímulo excitando um axônio, segundo o modelo de Fitzhugh-Nagumo. Essa simulação foi feita para <math>\alpha = 0.1</math> e um estímulo inicial de largura <math>\ell = 10</math> e intensidade <math>I=0.2</math> (acima do limiar). Note que quando o estímulo está acima do limiar, ele gera um PA que se propaga ao longo do axônio.]]

Nas animações das figuras 12 e 13 vemos a evolução espacial do potencial elétrico com o decorrer do tempo.
Na figura 12, o estímulo inicial está abaixo do limiar de voltagem, desta forma ele não gera um PA e a célula volta rapidamente para o seu PR. Na figura 13 o estímulo inicial está acima do limiar, logo um PA é gerado e o pulso elétrico é propagado ao longo do axônio.

===Modelo FitzHung-Nagumo 2D===

As simulação 2D tinham por objetivo simular a propagação de um estimulo inicial <math>I = 0.2</math> de <math>\ell = 5\times 50</math> de espessura. Para se aproximar da simetria do problema(axônio como um cilindro por onde o PA se propaga) aplicamos condições periódicas na vertical, e mantivemos as bordas nulas na horizontal.

Os parâmetros adotados para <math>\alpha, \epsilon, \gamma, D ,</math> e <math> \Delta x </math> foram os mesmos da secção 5.2. O <math> \Delta t</math> foi tomado como <math> 0.1 </math> para uma maior velocidade de processamento.

[[Arquivo:PA_Sob_Membrana2D.png|800px|thumb|center|'''Figura 14''' — Imagens de simulação de PA propagando sobre membrana 2D segundo modelo de FitzHung-Nagumo]]

Como podemos ver o PA se propaga assumindo um front de alta voltagem(despolarização e repolarização) estreito, e logo após uma região mais esparsa(de roxo mais escuro) indicando a repolarização da membrana da célula.

==Discussão==

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 <ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>

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

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 <math>\alpha > 0</math> <ref name=QUILLFELDT> https://www.ufrgs.br/mnemoforos/arquivos/potenciais2005.pdf Jorge A. Quillfeldt,"ORIGEM DOS POTENCIAIS ELÉTRICOS DAS CÉLULAS NERVOSAS" </ref>,que como mencionado na seçã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 <math>\alpha < 0</math>, mas embora essa abordagem leve a uma gama maior de resultados para estímulo contínuo
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>
<ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and Biology Society, Aug 2014, Chicago, United States. ffhal-00998828</ref> o sentido fisiológico é prejudicado, principalmente pelo fato de não possuir um limiar
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>.

==Programas==

===Cranck Nicolson e FTCS Modelo FitzHung-Nagumo 1D===
<source lang="python">

import numpy as np
import matplotlib
from matplotlib import pyplot as plt
from matplotlib.animation import ArtistAnimation

#Definindo Polinômio de Nagumo
def Pol(u):
a = 0.1
return u*(1 - u)*(u - a)

#Declarando Constantes
L = 100
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)

#Definindo função que calcula a equação de recuperação 1D
def FitzHung_FTCS(v, w):
D = 0.1
dx = 1
r = (D*dt)/(dx**2)
L = len(v)
aux = w[:]
for i in range(1, L-1):
aux[i] = w[i] + (r/(D*2))*0.005*(v[i] - 0.0001*w[i])

for j in range(1, L-1):
w[j] = aux[j]

return (w)

#Criando matriz tridiagonal para o método de Thomas
for i in range(L):
for j in range(L):
if i == j:
A[i,j] = 1 - 2*r
if j == i+1:
A[i,j] = r
if j == i-1:
A[i,j] = r

#Definindo função que calcula Equação de Nagumo 1D por método de Thomas
def Nagumo_CN(d,w,A):
L = len(d)
d = np.array(d)
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)

aux = np.dot(A, d)
aux[0] = d[0]
aux[L-1] = d[L-1]
for i in range(L):
d[i] = aux[i] + (r*Pol(aux[i]))/(D*2) - (r*w[i])/(2*D) + (r/(2*D))*I[i]

a = (-1*r)*(np.ones(L, dtype = np.float32))

b = (2*r + 1)*(np.ones(L, dtype = np.float32))

c = (-1*r)*(np.ones(L, dtype = np.float32))
c[L-1] = 0
g = [x for x in d]

c[0] = c[0]/b[0]
for i in range(1,L-1) :
c[i] = c[i]/(b[i] - (c[i-1]*a[i]))

for i in range(1,L-1) :
d[i] = (d[i] - (d[i-1]*a[i]))/(b[i] - (c[i-1]*a[i]))

g[L-1] = d[L-1]
for i in range(L-2,0,-1) :
g[i] = (d[i] - c[i]*g[i+1])

return (g)

#calculando o modelo FN
f = np.zeros(L, dtype = np.float32)
w = np.zeros(L, dtype = np.float32)
f[1:11] = 0.2

h = np.zeros(30000, dtype = np.float32)

T = np.arange(30000, dtype = np.float32)

x = np.arange(L, dtype = np.float32)
t = 0
h[0] = f[15]
for m in range(1, 30000):

aux = f[:]
f = Nagumo_CN(f, w, A)
w = FitzHung_FTCS(aux , w)
h[m] = f[15]
t = m*dt

#plot das figuras 10 e 11

ones = np.ones(30000)
zeros = np.zeros(30000)
plt.title('Modelo FN, $\\upsilon \\times t$ ($\ell = 10$, $I = 0.2$)')
plt.xlabel('t')
plt.ylabel('$\\upsilon$')
plt.plot(T, h)
plt.plot(T, ones,'r--', linewidth=0.8)
plt.plot(T, zeros, 'k--', linewidth = 0.8)
plt.savefig('nomedoarquivo.png')

#plot das animações

f = np.zeros(L, dtype = np.float32)
w = np.zeros(L, dtype = np.float32)
f[1:11] = 0.2
x = np.arange(L, dtype = np.float32)
t = 0
n = 1

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
frames = []
plt.title('Model FN, $\\ell$ = 10, I = 0.2')
plt.ylabel('$\\upsilon$')
plt.xlabel('x')
for m in range(1, 10000):

if (t == n):
curve = ax.plot(x, f, 'c')
frames.append(curve)
n = n + 1
aux = f[:]
f = Nagumo_CN(f, w, A)
w = FitzHung_FTCS(aux , w)

t = m*dt

animacao = ArtistAnimation(fig, frames, interval=50, blit=True)
animacao.save('nomedoarquivo.mp4')

</source>

===FTCS Modelo FitzHung-Nagumo 2D===
<source lang="python">

import numpy as np
import matplotlib
from matplotlib import pyplot as plt
import seaborn as sb

#Definindo Polinômio de Nagumo
def Pol(u):
a = 0.1
return u*(1 - u)*(u - a)

#Declarando Constantes
L = 100
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)

#Definindo função que calcula a equação de recuperação 2D
def FitzHung2D_FTCS(v, w):
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)
L = len(v)

aux = w[:]
for i in range(1, L-1):
for j in range(1, L-1):
aux[i,j] = w[i,j] + dt*0.005*(v[i,j] - 0.0001*w[i,j])

for i in range(1, L-1):
for j in range(1, L-1):
w[i,j] = aux[i,j]

return (w)

#Definindo função que calcula a equação de Nagumo em 2D como condições periódicas em y
def Nagumo2D_FTCS(v, w):
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)
L = len(v)

aux = v[:,:]
for j in range(1,L-1):
aux[0,j] = v[0,j] + r*(v[ L-1 ,j] + v[1,j] + v[0,j-1] + v[0,j+1] -4*v[0,j]) + dt*Pol(v[0,j]) - dt*w[0,j]

aux[0,L-1] = v[0,L-1] + r*(v[ L-1 ,L-1] + v[1,L-1] + v[0,L-2] + 0 -4*v[0,L-1]) + dt*Pol(v[0,L-1]) - dt*w[0,L-1]

for i in range(1, L-1):
for j in range(1, L-1):
aux[i,j] = v[i,j] + r*(v[i-1,j] + v[i+1,j] + v[i,j-1] + v[i,j+1] -4*v[i,j]) + dt*Pol(v[i,j]) - dt*w[i,j]

for j in range(1,L-1):
aux[L-1,j] = v[L-1,j] + r*(v[L-2,j] + v[0,j] + v[L-1,j-1] + v[L-1,j+1] -4*v[L-1,j]) + dt*Pol(v[L-1,j]) - dt*w[L-1,j]

aux[L-1,L-1] = v[L-1,L-1] + r*(v[L-2,L-1] + v[0,L-1] + v[L-1,L-2] + 0 -4*v[L-1,L-1]) + dt*Pol(v[L-1,L-1]) - dt*w[L-1,L-1]

for i in range(1, L):
for j in range(1, L-1):
v[i,j] = aux[i,j]

return (v)

#Função criada para evitar o uso de estruturas booleanas dentro do loop, que torna a simulação mais demorada
def Para_Evitar_Ifs(g):
intervalo = 20
T = g[0]
v = g[1]
w = g[2]
for t in range(T,T + intervalo):
aux = v[:,:]
v = Nagumo2D_FTCS(v, w)
w = FitzHung2D_FTCS(aux , w)

return ([T+intervalo,v,w])

#Criando os espaço e o grid de fotos
x = np.arange(L, dtype = np.float32)
y = np.arange(L, dtype = np.float32)
T = 0
dl = 4
p1 = (dl/2) - 1
p2 = (dl/2) + 1
tamanho = 2*L/dl

x,y = np.meshgrid(x,y)

f, (ax1,ax2,ax3,ax4) = plt.subplots(1,4,sharey=True, figsize=(24,4))
f.suptitle('Pulso Lateral de Tamanho '+str(int(tamanho))+'')

#Declarando condições iniciais e condições de contorno não periódicas
v = np.zeros((L,L))
v[int(p1*tamanho/2)-6:int(p2*tamanho/2),1:5] = 0.2

w = np.zeros((L,L))

#Plotando os 4 primeiros mapas de calor
g = [T,v,w]
graph1 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax1)
graph1.set_title('t='+str(g[0])+'')
graph1.invert_yaxis()

g = Para_Evitar_Ifs(g)
graph2 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax2)
graph2.set_title('t='+str(g[0])+'')
graph2.invert_yaxis()

g = Para_Evitar_Ifs(g)
graph3 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax3)
graph3.set_title('t='+str(g[0])+'')
graph3.invert_yaxis()

g = Para_Evitar_Ifs(g)
graph4 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax4)
graph4.set_title('t='+str(g[0])+'')
graph4.invert_yaxis()

plt.savefig('Lateral_1.png')
plt.show()

</source>

==Referências==
<references/>

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

2021-04-05T14:23:41Z

Nataliaferrazzo:

'''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)<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| A simple spiking neuron model: sodium channels alone </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 Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </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>
<ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </ref>
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>

<center><math>
\begin{cases}
\frac{\partial v}{\partial t} = D\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)
</math></center>

Onde <math>0<\alpha<1</math> é o limiar de potencial, <math>D</math> é uma constante de difusão, <math>\epsilon</math> e <math>\gamma</math> são parâmetros positivos a serem ajustados —relacionados à velocidade de atuação da variável de recuperação— e <math>I</math> é a magnitude do estímulo externo <ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </ref>
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>

== Método Crank-Nicolson ==

Uma equação diferencial parcial da forma
<center><math>Au_{xx} + Bu_{xy} + Cu_{yy} = f(x, y, u_x, u_y) \qquad (3.1)</math></center>

onde A, B e C são constantes, é chamada de equação quasilinear.
Um método implícito no tempo e estável numericamente de solucionar este tipo de equação foi proposto na metade do século XX por John Crank e Phyllis Nicolson.
O método de Crank-Nicolson é de segunda ordem no tempo e no espaço e é baseado em diferenças centradas no espaço e regra trapezoidal no tempo.
Como provaremos mais adiante, este método é incondicionalmente estável.

A equação diferencial finita de Nagumo é dada pela segunda derivada espacial obtida a partir de uma combinação da derivada temporal no passo <math> n</math> e <math> n+1</math>:

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} = (\theta) \frac{u_{i+1}^n -2u_i^{n} + u_{i+1}^n}{\Delta x^2} + (1-\theta) \frac{u_{i+1}^{n+1} -2u_i^{n+1} + u_{i-1}^{n+1}}{\Delta x^2} \qquad(3.2)</math></center>

A partir da equação (3.2), podemos obter os diferentes métodos de diferenciação da equação de Nagumo ajustando o valor de <math>\theta</math>.
*<math>\theta = 1</math> é o método explícito.
*<math>\theta = 0</math> é o método implícito.
*<math>\theta = \frac{1}{2}</math> é o método Crank-Nicolson.

Substituindo o valor de <math>\theta</math> na equação (3.2), temos:

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} = \frac{u_{i+1} ^n -2u_i ^{n} + u_{i+1} ^n}{\Delta x^2} + u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.3)</math></center>

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} = \frac{u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1}}{\Delta x^2} + u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha) \qquad(3.4)</math></center>

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} =\frac{u_{i+1} ^n -2u_i ^{n} + u_{i+1} ^n}{2\Delta x^2} + \frac{u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1}}{2\Delta x^2} + u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.5)</math></center>

Para nosso problema, usaremos o método pseudo Crank-Nicolson, devido ao termo não-linear presente em (3.5).

=== Análise de Estabilidade de Von Neumann ===
Baseada nas deduções de G. Garcia <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>.

Ao analisar a estabilidade do método, precisamos verificar se o erro gerado no procedimento não aumenta quando testamos os extremos da solução.
Para realizar esta verificação, vamos utilizar o método de Von Neumann.
A primeira coisa a se fazer é linearizar o termo não linear da equação (3.5):

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} =\frac{u_{i+1} ^n -2u_i ^{n} + u_{i+1} ^n}{2\Delta x^2} + \frac{u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1}}{2\Delta x^2} -\alpha u_{i} ^{n}\qquad(3.6)</math></center>

Se substituírmos <math>u_{i} ^{n} = \Psi (t)e^{j\gamma \Delta x}</math></div> no tempo <math>t</math>, onde <math>j = \sqrt{(-1)}</math> e <math>\gamma > 0</math> o método se mostra estável pelo seguinte:

<center><math>\bigg| \frac{\Psi (t + \Delta t)}{\Psi (t)} \bigg| \le 1 </math></center>

Substituindo <math>u_{i} ^{n} = \Psi (t)e^{j\gamma x} </math> na equação (3.6), temos:

<center><math>\begin{matrix} \frac{\Psi (t + \Delta t) e^{j\gamma x}}{\Delta t}& =&\frac{\Psi (t) e^{j\gamma (x + \Delta x)} -2 \Psi (t) e^{j\gamma x} + \Psi (t) e^{j\gamma (x - \Delta x)}}{2\Delta x^2} \\ \ &+ &\frac{\Psi (t + \Delta t) e^{j\gamma (x + \Delta x)} -2\Psi (t + \Delta t)e^{j\gamma x} + \Psi (t + \Delta t) e^{j\gamma (x + \Delta x)}}{2\Delta x^2} -\alpha \Psi (t) e^{j\gamma x}\end{matrix}\qquad(3.7)</math></center>

Substituindo <math>r = {\Delta t}/{\Delta x^2}</math> em (3.7), temos:

<center><math>\begin{matrix} \Psi (t + \Delta t) e^{j\gamma x} &=& r \bigg(\frac{\Psi (t) e^{j\gamma (x + \Delta x)} -2\Psi (t) e^{j\gamma x} + \Psi (t) e^{j\gamma (x - \Delta x)}}{2}\bigg) \\ \ &
+& r\bigg( \frac{\Psi (t + \Delta t) e^{j\gamma (x + \Delta x)} -2\Psi (t + \Delta t)e^{j\gamma x} + \Psi (t + \Delta t) e^{j\gamma (x + \Delta x)}}{2}\bigg) -\Delta t \alpha \Psi (t) e^{j\gamma x}\end{matrix}\qquad(3.8)</math></center>

Podemos notar que há <math>e^{j\gamma x} </math> em todos os termos. Simplificando-os e movendo todos os termos com <math>\Psi (t + \Delta t) </math> para o lado esquerdo e todos com <math>\Psi (t) </math> para a direita, a equação (3.8) toma a forma:

<center><math>\Psi (t + \Delta t) - r \bigg(\frac{\Psi (t + \Delta t)e^{j\gamma \Delta x} -2\Psi (t + \Delta t) + \Psi (t + \Delta t)e^{-j\gamma \Delta x}}{2}\bigg) = r \bigg( \frac{\Psi (t) e^{j\gamma \Delta x} -2\Psi (t) + \Psi (t) e^{-j\gamma \Delta x}}{2}\bigg) + \Psi (t) -\Delta t \alpha \Psi (t)\qquad(3.9)</math></center>

Fatorando os termos, ficamos com:

<center><math>\Psi (t + \Delta t) \bigg( 1 - r \bigg( \frac{e^{j\gamma \Delta x} -2 + e^{-j\gamma \Delta x}}{2}\bigg) \bigg) = \Psi (t) \bigg( r \bigg( \frac{e^{j\gamma \Delta x} - 2 + e^{-j\gamma \Delta x}}{2} \bigg) + 1 -\Delta t \alpha \bigg) \qquad(3.10)</math></center>

Utilizando a identidade trigonométrica <math> cos(x) = \frac{e^{jx} + e^{-jx}}{2}</math>, ficamos com:

<center><math>\Psi (t + \Delta t) (1 - r (cos(\gamma \Delta x) -1)) = \Psi (t) (r (cos(j\gamma \Delta x) -1 ) +1 -\Delta t \alpha)\qquad(3.11)</math></center>

Aplicando <math>1 - cos(2u) = 2 sin^2 (u)</math>, temos:

<center><math>\Psi (t + \Delta t) \bigg(1 + 2r sin^2 \bigg( \frac{\gamma \Delta x}{2}\bigg) \bigg) = \Psi (t) \bigg(1 - 2r sin^2 \bigg( \frac{\gamma \Delta x}{2}\bigg) -\Delta t \alpha\bigg)\qquad(3.12)</math></center>

Pode-se demonstrar que, a partir da equação (3.12), tem-se:

<center><math>\frac{\Psi (t + \Delta t) e^{j\gamma x}}{\Delta t} =\frac{1 - 2r sin^2( \frac{\gamma \Delta x}{2}) -\Delta t \alpha}{1 + 2r sin^2( \frac{\gamma \Delta x}{2})}\qquad(3.13)</math></center>

Como:

<center><math>\bigg| \frac{1 - 2r sin^2( \frac{\gamma \Delta x}{2}) -\Delta t \alpha}{1 + 2r sin^2( \frac{\gamma \Delta x}{2})}\bigg| \le 1 \qquad(3.13)</math></center>

prova-se que o método Crank-Nicolson é estável para <math>\Delta t </math> suficientemente pequenos.

=== Implementação para a equação de FitzHung-Nagumo ===

Para a implementação na equação de FitzHung-Nagumo, as seguintes derivadas parciais são utilizadas:

<center><math>u(x_{i} , t_{n}) = u_{i}^n \qquad(3.14)</math></center>

<center><math>\frac{\partial u}{\partial t} (x_{i} , t_{n})= \frac{u_{i}^{n+1}- u_{i}^{n}}{\Delta t}\qquad(3.15)</math></center>

<center><math>\frac{\partial^2 u}{\partial^2 x} (x_{i} , t_{n})= \frac{1}{2}\bigg(\frac{u_{i+1}^{n+1}- 2u_{i}^{n+1}+u_{i-1}^{n+1}}{\Delta x^2} + \frac{u_{i+1}^{n}- 2u_{i}^{n}+u_{i-1}^{n}}{\Delta x^2}\bigg)\qquad(3.16)</math></center>

Neste ponto, utilizamos uma aproximação para a derivada temporal e a derivada espacial de segunda ordem e substituímos as equações (3.14), (3.15) e (3.16) na equação (2).

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} =\frac{1}{2\Delta x^2}(u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1})+ \frac{1}{2\Delta x^2}(u_{i+1} ^{n} -2u_i ^{n} + u_{i-1} ^{n})
u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.17)</math></center>

Multiplicando por <math>\Delta t</math>, temos:

<center><math>u_i ^{n+1} - u_i ^{n} =\frac{\Delta t}{2\Delta x^2}(u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1})+ \frac{\Delta t}{2\Delta x^2}(u_{i+1} ^{n} -2u_i ^{n} + u_{i-1} ^{n})
\Delta t u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.18)</math></center>

Fazendo a substituição:

<center><math> r = \frac{2\Delta t}{\Delta x^2}</math></center>

e rearranjando, temos:

<center><math>-ru_{i+1}^{n+1} + (1+2r)u_{i}^{n+1} - ru_{i-1}^{n+1} = r(u_{i+1}^{n} + u_{i-1}^{n}) + (1-2r)u_{i}^{n} + \Delta tu_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.19)</math></center>

Podemos, então representar a equação (3.19) em forma matricial:

<center><math>\mathbf{AU} = \mathbf{F}</math></center>

Onde o lado esquerdo da equação é a matriz tridimensional a seguir
<center><math>
A =
\begin{bmatrix}
1+2r & -r & 0 & \cdots & 0 \\
-r & 1+2r & -r & \cdots & 0 \\
0 & -r & 1+2r & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & -r \\
0 & 0 & 0 & -r & 1+2r \\
\end{bmatrix}
</math></center>

e o lado direito:
<center><math>
F =
\begin{bmatrix}
0 \\
\vdots \\
r(u_{n+1 ^n} + u_{i-1 ^n}) + (1-2r)u_i ^n + \Delta t u_i ^n(1-u_i ^n)(u_i ^n - \alpha) \\
\vdots \\
1\\
\end{bmatrix}
</math></center>

Com isso, podemos encontrar as soluções da equação resolvendo:

<center><math>\mathbf{U} = \mathbf{F/A}</math></center>

===FTCS Explícito===
No método explícito, a discretização é

<math>
\begin{cases}
\frac{\partial f}{\partial t} \rightarrow \frac{f(x, t+\Delta t) - f(x,t)}{\Delta t}\\
\frac{\partial^2 f}{\partial x^2} \rightarrow \frac{f(x, t+\Delta t) - 2f(x,t) + f(x-\Delta t}{(\Delta x)^2}
\end{cases}
</math>

Substituindo na EDU e utilizando a notação de índices (<math>n</math> é índice temporal e <math>j</math> é o índice temporal), temos

<math>
\frac{f_j^{n+1} - f_j^n}{\Delta t} = D\frac{f_{j+1}^n - 2f_j^n + f_{j-1}^n}{(\Delta x)^2}
</math>

Agrupando <math>\Delta t,\ \Delta x\ e\ D</math> em uma constante <math>k</math> e deixando o passo temporal futuro em função do passo anterior explicitamente, temos

<math>
f_j^{n+1} = f_j^n+k(f_{j+1}^n - 2f_j^n + f_{j-1}^n).
</math>

Utilizando esta equação é possível calcular o valor de <math>f</math> iterativamente para todo intervalo de tempo.

===FTCS Implícito===
No método implícito, também chamado de BTCS (para trás no tempo centrado no espaço, do inglês: Backwards Time Centered Space), pois o passo temporal é dado para trás, assim a derivada temporal é

<math>
\begin{cases}
\frac{\partial f}{\partial t} \rightarrow \frac{f_j^{n} - f_j^{n-1}}{\Delta t}\\
\end{cases}
</math>

Resultando na seguinte EDU

<math>
f_j^{n} = f_j^{n-1}+k(f_{j+1}^n - 2f_j^n + f_{j-1}^n).
</math>

Reescrevendo os índices <math>n</math> temos

<math>
f_j^{n+1} = f_j^{n}+k(f_{j+1}^{n+1} - 2f_j^{n+1} + f_{j-1}^{n+1})
</math>

Onde <math>f_j^{n+1}</math> está escrita implicitamente, e pode ser calculada utilizando o Algoritmo de Thomas [https://moodle.ufrgs.br/pluginfile.php/3653465/mod_resource/content/1/thomas.pdf].

===Equação de Recuperação 1D===

No modelo de Fitzhugh-Nagumo, utilizou-se o mesmo algoritmo de Crank-Nicolson utilizado na Equação de Nagumo, na seção 3.1. Entretanto, ainda é necessário resolver a equação diferencial de <math>w</math>, dada por

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

Para isso utilizou-se o método FTCS explícito. Discretizando a derivada, temos

<math>
\frac{\partial w}{\partial t} \rightarrow \frac{w(x,t + \Delta t)-w(x,t)}{\Delta t}
</math>

Utilizando a notação de índices e aplicando na equação da recuperação, temos

<math>
\frac{w_j^{n+1}-w_j^n}{\Delta t} = \epsilon(v_j^n - \gamma w_j^n)
</math>

<math>
w_j^{n+1}-w_j^n = \Delta t \epsilon(v_j^n - \gamma w_j^n).
</math>

Assim, uma vez conhecido <math>v_j^n</math>, é possível calcular <math>w</math> para qualquer intervalo de tempo.

Vamos verificar as condições de estabilidade do método utilizando o os modos de Fourier:
*substituindo <math>w_j^n=A^ne^{i(jq\Delta x)},</math> onde <math>i</math> é a unidade imaginária;
*desconsiderando <math>v_j^n</math>,
temos

<math>
A^{n+1}e^{i(jq \Delta x)} = -\Delta t \epsilon \gamma e^{i(lq \delta x)}
</math>

Dividindo ambos os lados por <math>A^ne^{i(jq\Delta x)},\ </math> temos

<math>
\frac{A^{n+1}}{A^n} = -\Delta t \epsilon \gamma
</math>

Para que o método seja estável, temos que <math>\Bigg|\frac{A^{n+1}}{A^n}\Bigg| < 1,\ </math> logo

<math>
\Bigg|\frac{A^{n+1}}{A^n}\Bigg| =\Bigg|-\Delta t \epsilon \gamma \Bigg|< 1.
</math>

Assim, temos que a condição de estabilidade sobre o método é tal que

<math>
\frac{-1}{\epsilon \gamma}<\Delta t < \frac{1}{\epsilon \gamma}.
</math>

===Modelo FitzHung-Nagumo 2D===
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 de Nagumo em 2D. Assumindo <math>\Delta x = \Delta y</math>, temos

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

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

<math>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} - 4 v_{l,j}^{n}\right] + \Delta t f_{N l,j}^{n} - \Delta t w_{l,j}^{n}</math>

onde os índices <math>l</math> e <math>j</math> se referem às coordenadas <math>x</math> e <math>y</math>, respectivamente, e o índice <math>n</math> refere-se ao tempo.

Vamos verificar as condições de estabilidade do método utilizando os modos de Fourier:
*substituindo <math>v_{l,j}^{n} = A^{n}e^{i(xq_x +yq_y)} = A^{n}e^{i\beta}</math>
*usando <math>r = \frac{D \Delta t}{(\Delta x)^2}</math>
*por simplicidade, como fizemos na seção 3.1, vamos linearizar <math>f_{N l,j}^{n} \rightarrow -\alpha v_{l,j}^{n} =-\alpha A^{n}e^{i\beta}</math>
*por fim, desconsiderando <math>\Delta t w_{l,j}^{n}</math>, temos

<math>A^{n+1}e^{i\beta} = A^{n}e^{i\beta} + rA^{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)} - 4 e^{i\beta}\right] -\Delta t \alpha A^{n}e^{i\beta}</math>

<math>A^{n+1}e^{i\beta} = A^{n}e^{i\beta} + rA^{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}</math>

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

<math>\frac{A^{n+1}}{A^{n}} = 1 + r\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</math>

Sabendo as seguintes identidades trigonométricas:

*<math>e^{-i\theta} + e^{i\theta} = 2\cos(\theta)</math>
*<math>\cos(2\theta) = 1 - 2\sin^2(\theta)</math>

e aplicando à equação, obtemos

<math>\frac{A^{n+1}}{A^{n}} = 1 + 2r\left[ \cos(\Delta xq_x) + \cos(\Delta yq_y) - 2\right] -\Delta t \alpha</math>

<math>\frac{A^{n+1}}{A^{n}} = 1 + 2r\left[ (1 - 2\sin^2(\Delta xq_x/2)) + (1 - 2\sin^2(\Delta yq_y/2)) - 2\right] -\Delta t \alpha</math>

<math>\frac{A^{n+1}}{A^{n}} = 1 + 4r\left[ \sin^2(\Delta xq_x/2) + \sin^2(\Delta yq_y/2)\right] -\Delta t \alpha.</math>

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

<math>\left|1 -\Delta t \alpha + 4r\left[ \sin^2(\Delta xq_x/2) + \sin^2(\Delta yq_y/2)\right]\right| < 1</math>.

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

<math>\frac{\alpha \Delta t}{8} < r < \frac{2 + \alpha \Delta t}{8}</math>.

==Resultados==
===Equação de Nagumo 1D===

Para as simulações numéricas foram utilizados os seguintes parâmetros:

* <math>\alpha = 0,1</math>
* <math>\Delta t = 0,01\qquad para 0< t < 20000</math>
* <math>\Delta x = 0,01\qquad para 0< x < 200</math>
* <math>\Delta t = 0,01\qquad para 0<\Delta t < 20000</math>

Com isso, dividimos as simulações em dois casos.

Caso I: As simulações foram geradas a partir da definição de uma largura de estímulo estreita <math>l= 20</math>.
Com isso, variamos a intensidade do estímulo. Nas seguintes figuras, podemos ver as simulações para uma intensidade menor do que o limiar, igual ao limiar e maior que o limiar, respectivamente, <math>I < \alpha</math>, <math>I = \alpha</math> e <math>I > \alpha</math>.

[[Arquivo:N_fino_abaixolimiar.gif|800px|thumb|center|'''Figura 4''' — estímulo <math>\ell = 20</math> <math>I = 0.05</math>]]

Como pode ser visto na Figura 4, quando <math>I=0,05</math>, ou seja, menor que o limiar, o sinal rapidamente se esvai.

Para uma intensidade igual ao limiar, o sinal também se extingue, porém com um intervalo de tempo maior do que anteriormente.
[[Arquivo:N_fino_iguallimiar.gif|800px|thumb|center|'''Figura 6''' — estímulo <math>\ell = 20</math> <math>I = 0.1</math>]]

E, finalmente, quando a intensidade é maior que o limiar, não importa quão maior seja, o sinal é passado ao máximo da função.
[[Arquivo:N_fino_acimalimiar.gif|800px|thumb|center|'''Figura 5''' — estímulo <math>\ell = 20</math> <math>I = 0.2</math>]]

Caso II: As seguintes simulações foram geradas a partir da definição de uma largura de estímulo larga <math>l= 100</math>.
Com isso, novamente a intensidade do estímulo é variada.

[[Arquivo:N_largo_abaixolimiar.gif|800px|thumb|center|'''Figura 7''' — estímulo <math>\ell = 100</math> <math>I = 0.05</math>]]
Para uma intensidade menor que o limiar, novamente o sinal se esvai.

[[Arquivo:N_largo_iguallimiar.gif|800px|thumb|center|'''Figura 9''' — estímulo <math>\ell = 100</math> <math>I = 0.1</math>]]
O mesmo acontece quando a intensidade é igual ao limiar.

Já quando a intensidade é maior que o limiar, o comportamento “tudo ou nada” novamente se mostra presente.
[[Arquivo:N_largo_acimalimiar.gif|800px|thumb|center|'''Figura 8''' — estímulo <math>\ell = 100</math> <math>I = 0.2</math>]]

Podemos notar que quando a largura é maior, o sinal se esvai mais lentamente se comparado a uma largura de sinal mais estreita.

===Modelo FitzHung-Nagumo 1D===
Para a implementação da simulação do modelo FN em 1D utilizou-se <math>\Delta t = 0.01</math> (com <math>0<t<30000</math>), <math>\Delta x = 1</math> (com <math>0<x<120</math>) e os demais parâmetros foram idênticos àqueles encontrados na ref.[1] <ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </ref>. Eles são

<math>
\begin{cases}
\alpha = 0.1\\
\epsilon = 0.005\\
\gamma = 0.0001\\
D = 0.1\\
\end{cases}
</math>

Foram aplicadas condições de contorno nulas nas bordas do intervalo de <math>x</math>.

[[Arquivo:FN_vXt.png|800px|thumb|center|'''Figura 10''' — Evolução temporal do potencial elétrico em um ponto da membrana do axônio quando excitado por um estímulo acima do limiar de potencial.]]

Na figura 10 vemos a evolução temporal do potencial elétrico em um ponto da membrana do axônio (<math>x=15</math>) quando excitado acima do limiar de potencial. Note que a curva do PA é muito semelhante àquela da na Figura 2, apresentando vários dos aspectos característicos, como a assimetria entre os períodos de polarização (ascensão rápida) e repolarização (decaimento mais lento) e um período refratário, onde o potencial lentamente se recupera em direção do PR. Podemos assim entender a escolha dos valores dos parâmetros <math>\epsilon</math> e <math>\gamma</math>, que controlam a velocidade de crescimento da variável de recuperação (assim como se pode ver na equação (3)). Experimentalmente, sabemos que o período refratário é muito mais longo do que o período de despolarização
<ref name=QUILLFELDT> https://www.ufrgs.br/mnemoforos/arquivos/potenciais2005.pdf Jorge A. Quillfeldt,"ORIGEM DOS POTENCIAIS ELÉTRICOS DAS CÉLULAS NERVOSAS" </ref>
, logo precisamos que a intensidade da variável <math>w</math> aumente muito lentamente. Para isso, escolhe-se valores pequenos de <math>\epsilon</math> e <math>\gamma</math>, pois caso contrário o potencial da membrana celular não teria tempo de atingir o pico do PA.

A velocidade de crescimento da variável de recuperação também explica a ausência do "tudo ou nada" no PA da figura 10. Na figura 11, vê-se a evolução temporal de um PA com <math>\epsilon = 0.001</math> (5 vezes menor do que aquele da figura 10).

[[Arquivo:FN_vXt_Epequeno.png|800px|thumb|center|'''Figura 11''' — Evolução temporal do potencial elétrico em um ponto da membrana do axônio para <math>\epsilon = 0.001</math>.]]

Na figura 11, como a variável de recuperação cresce muito mais lentamente, além de demorar mais tempo para a membrana celular sair do período refratário, vemos que o pico do PA de ação é muito mais próximo do valor máximo. Isso ocorre porque a variável de recuperação cresce mais lentamente, possibilitando que a variável rápida atinja um valor maior antes de começar a diminuir devido à ação de <math>w</math>.

[[Arquivo:FN_abaixolim_apos.gif|800px|thumb|center|'''Figura 12''' — Simulação de um estímulo excitando um axônio, segundo o modelo de Fitzhugh-Nagumo. Essa simulação foi feita para <math>\alpha = 0.1</math> e um estímulo inicial de largura <math>\ell = 10</math> e intensidade <math>I=0.05</math> (abaixo do limiar). Note que o estímulo não desencadeia um PA e logo desaparece.]]

[[Arquivo:FN_acimalim_apos.gif|800px|thumb|center|'''Figura 13''' — Simulação de um estímulo excitando um axônio, segundo o modelo de Fitzhugh-Nagumo. Essa simulação foi feita para <math>\alpha = 0.1</math> e um estímulo inicial de largura <math>\ell = 10</math> e intensidade <math>I=0.2</math> (acima do limiar). Note que quando o estímulo está acima do limiar, ele gera um PA que se propaga ao longo do axônio.]]

Nas animações das figuras 12 e 13 vemos a evolução espacial do potencial elétrico com o decorrer do tempo.
Na figura 12, o estímulo inicial está abaixo do limiar de voltagem, desta forma ele não gera um PA e a célula volta rapidamente para o seu PR. Na figura 13 o estímulo inicial está acima do limiar, logo um PA é gerado e o pulso elétrico é propagado ao longo do axônio.

===Modelo FitzHung-Nagumo 2D===

As simulação 2D tinham por objetivo simular a propagação de um estimulo inicial <math>I = 0.2</math> de <math>\ell = 5\times 50</math> de espessura. Para se aproximar da simetria do problema(axônio como um cilindro por onde o PA se propaga) aplicamos condições periódicas na vertical, e mantivemos as bordas nulas na horizontal.

Os parâmetros adotados para <math>\alpha, \epsilon, \gamma, D ,</math> e <math> \Delta x </math> foram os mesmos da secção 5.2. O <math> \Delta t</math> foi tomado como <math> 0.1 </math> para uma maior velocidade de processamento.

[[Arquivo:PA_Sob_Membrana2D.png|800px|thumb|center|'''Figura 14''' — Imagens de simulação de PA propagando sobre membrana 2D segundo modelo de FitzHung-Nagumo]]

Como podemos ver o PA se propaga assumindo um front de alta voltagem(despolarização e repolarização) estreito, e logo após uma região mais esparsa(de roxo mais escuro) indicando a repolarização da membrana da célula.

==Discussão==

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 <ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>

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

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 <math>\alpha > 0</math> <ref name=QUILLFELDT> https://www.ufrgs.br/mnemoforos/arquivos/potenciais2005.pdf Jorge A. Quillfeldt,"ORIGEM DOS POTENCIAIS ELÉTRICOS DAS CÉLULAS NERVOSAS" </ref>,que como mencionado na seçã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 <math>\alpha < 0</math>, mas embora essa abordagem leve a uma gama maior de resultados para estímulo contínuo
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>
<ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and Biology Society, Aug 2014, Chicago, United States. ffhal-00998828</ref> o sentido fisiológico é prejudicado, principalmente pelo fato de não possuir um limiar
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>.

==Programas==

===Cranck Nicolson e FTCS Modelo FitzHung-Nagumo 1D===
<source lang="python">

import numpy as np
import matplotlib
from matplotlib import pyplot as plt
from matplotlib.animation import ArtistAnimation

#Definindo Polinômio de Nagumo
def Pol(u):
a = 0.1
return u*(1 - u)*(u - a)

#Declarando Constantes
L = 100
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)

#Definindo função que calcula a equação de recuperação 1D
def FitzHung_FTCS(v, w):
D = 0.1
dx = 1
r = (D*dt)/(dx**2)
L = len(v)
aux = w[:]
for i in range(1, L-1):
aux[i] = w[i] + (r/(D*2))*0.005*(v[i] - 0.0001*w[i])

for j in range(1, L-1):
w[j] = aux[j]

return (w)

#Criando matriz tridiagonal para o método de Thomas
for i in range(L):
for j in range(L):
if i == j:
A[i,j] = 1 - 2*r
if j == i+1:
A[i,j] = r
if j == i-1:
A[i,j] = r

#Definindo função que calcula Equação de Nagumo 1D por método de Thomas
def Nagumo_CN(d,w,A):
L = len(d)
d = np.array(d)
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)

aux = np.dot(A, d)
aux[0] = d[0]
aux[L-1] = d[L-1]
for i in range(L):
d[i] = aux[i] + (r*Pol(aux[i]))/(D*2) - (r*w[i])/(2*D) + (r/(2*D))*I[i]

a = (-1*r)*(np.ones(L, dtype = np.float32))

b = (2*r + 1)*(np.ones(L, dtype = np.float32))

c = (-1*r)*(np.ones(L, dtype = np.float32))
c[L-1] = 0
g = [x for x in d]

c[0] = c[0]/b[0]
for i in range(1,L-1) :
c[i] = c[i]/(b[i] - (c[i-1]*a[i]))

for i in range(1,L-1) :
d[i] = (d[i] - (d[i-1]*a[i]))/(b[i] - (c[i-1]*a[i]))

g[L-1] = d[L-1]
for i in range(L-2,0,-1) :
g[i] = (d[i] - c[i]*g[i+1])

return (g)

#calculando o modelo FN
f = np.zeros(L, dtype = np.float32)
w = np.zeros(L, dtype = np.float32)
f[1:11] = 0.2

h = np.zeros(30000, dtype = np.float32)

T = np.arange(30000, dtype = np.float32)

x = np.arange(L, dtype = np.float32)
t = 0
h[0] = f[15]
for m in range(1, 30000):

aux = f[:]
f = Nagumo_CN(f, w, A)
w = FitzHung_FTCS(aux , w)
h[m] = f[15]
t = m*dt

#plot das figuras 10 e 11

ones = np.ones(30000)
zeros = np.zeros(30000)
plt.title('Modelo FN, $\\upsilon \\times t$ ($\ell = 10$, $I = 0.2$)')
plt.xlabel('t')
plt.ylabel('$\\upsilon$')
plt.plot(T, h)
plt.plot(T, ones,'r--', linewidth=0.8)
plt.plot(T, zeros, 'k--', linewidth = 0.8)
plt.savefig('nomedoarquivo.png')

#plot das animações

f = np.zeros(L, dtype = np.float32)
w = np.zeros(L, dtype = np.float32)
f[1:11] = 0.2
x = np.arange(L, dtype = np.float32)
t = 0
n = 1

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
frames = []
plt.title('Model FN, $\\ell$ = 10, I = 0.2')
plt.ylabel('$\\upsilon$')
plt.xlabel('x')
for m in range(1, 10000):

if (t == n):
curve = ax.plot(x, f, 'c')
frames.append(curve)
n = n + 1
aux = f[:]
f = Nagumo_CN(f, w, A)
w = FitzHung_FTCS(aux , w)

t = m*dt

animacao = ArtistAnimation(fig, frames, interval=50, blit=True)
animacao.save('nomedoarquivo.mp4')

</source>

===FTCS Modelo FitzHung-Nagumo 2D===
<source lang="python">

import numpy as np
import matplotlib
from matplotlib import pyplot as plt
import seaborn as sb

#Definindo Polinômio de Nagumo
def Pol(u):
a = 0.1
return u*(1 - u)*(u - a)

#Declarando Constantes
L = 100
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)

#Definindo função que calcula a equação de recuperação 2D
def FitzHung2D_FTCS(v, w):
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)
L = len(v)

aux = w[:]
for i in range(1, L-1):
for j in range(1, L-1):
aux[i,j] = w[i,j] + dt*0.005*(v[i,j] - 0.0001*w[i,j])

for i in range(1, L-1):
for j in range(1, L-1):
w[i,j] = aux[i,j]

return (w)

#Definindo função que calcula a equação de Nagumo em 2D como condições periódicas em y
def Nagumo2D_FTCS(v, w):
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)
L = len(v)

aux = v[:,:]
for j in range(1,L-1):
aux[0,j] = v[0,j] + r*(v[ L-1 ,j] + v[1,j] + v[0,j-1] + v[0,j+1] -4*v[0,j]) + dt*Pol(v[0,j]) - dt*w[0,j]

aux[0,L-1] = v[0,L-1] + r*(v[ L-1 ,L-1] + v[1,L-1] + v[0,L-2] + 0 -4*v[0,L-1]) + dt*Pol(v[0,L-1]) - dt*w[0,L-1]

for i in range(1, L-1):
for j in range(1, L-1):
aux[i,j] = v[i,j] + r*(v[i-1,j] + v[i+1,j] + v[i,j-1] + v[i,j+1] -4*v[i,j]) + dt*Pol(v[i,j]) - dt*w[i,j]

for j in range(1,L-1):
aux[L-1,j] = v[L-1,j] + r*(v[L-2,j] + v[0,j] + v[L-1,j-1] + v[L-1,j+1] -4*v[L-1,j]) + dt*Pol(v[L-1,j]) - dt*w[L-1,j]

aux[L-1,L-1] = v[L-1,L-1] + r*(v[L-2,L-1] + v[0,L-1] + v[L-1,L-2] + 0 -4*v[L-1,L-1]) + dt*Pol(v[L-1,L-1]) - dt*w[L-1,L-1]

for i in range(1, L):
for j in range(1, L-1):
v[i,j] = aux[i,j]

return (v)

#Função criada para evitar o uso de estruturas booleanas dentro do loop, que torna a simulação mais demorada
def Para_Evitar_Ifs(g):
intervalo = 20
T = g[0]
v = g[1]
w = g[2]
for t in range(T,T + intervalo):
aux = v[:,:]
v = Nagumo2D_FTCS(v, w)
w = FitzHung2D_FTCS(aux , w)

return ([T+intervalo,v,w])

#Criando os espaço e o grid de fotos
x = np.arange(L, dtype = np.float32)
y = np.arange(L, dtype = np.float32)
T = 0
dl = 4
p1 = (dl/2) - 1
p2 = (dl/2) + 1
tamanho = 2*L/dl

x,y = np.meshgrid(x,y)

f, (ax1,ax2,ax3,ax4) = plt.subplots(1,4,sharey=True, figsize=(24,4))
f.suptitle('Pulso Lateral de Tamanho '+str(int(tamanho))+'')

#Declarando condições iniciais e condições de contorno não periódicas
v = np.zeros((L,L))
v[int(p1*tamanho/2)-6:int(p2*tamanho/2),1:5] = 0.2

w = np.zeros((L,L))

#Plotando os 4 primeiros mapas de calor
g = [T,v,w]
graph1 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax1)
graph1.set_title('t='+str(g[0])+'')
graph1.invert_yaxis()

g = Para_Evitar_Ifs(g)
graph2 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax2)
graph2.set_title('t='+str(g[0])+'')
graph2.invert_yaxis()

g = Para_Evitar_Ifs(g)
graph3 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax3)
graph3.set_title('t='+str(g[0])+'')
graph3.invert_yaxis()

g = Para_Evitar_Ifs(g)
graph4 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax4)
graph4.set_title('t='+str(g[0])+'')
graph4.invert_yaxis()

plt.savefig('Lateral_1.png')
plt.show()

</source>

==Referências==
<references/>

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

2021-04-05T02:34:47Z

Nataliaferrazzo:

'''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)<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| A simple spiking neuron model: sodium channels alone </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 Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </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>
<ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </ref>
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>

<center><math>
\begin{cases}
\frac{\partial v}{\partial t} = D\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)
</math></center>

Onde <math>0<\alpha<1</math> é o limiar de potencial, <math>D</math> é uma constante de difusão, <math>\epsilon</math> e <math>\gamma</math> são parâmetros positivos a serem ajustados —relacionados à velocidade de atuação da variável de recuperação— e <math>I</math> é a magnitude do estímulo externo <ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </ref>
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>

== Método Crank-Nicolson ==

Uma equação diferencial parcial da forma
<center><math>Au_{xx} + Bu_{xy} + Cu_{yy} = f(x, y, u_x, u_y) \qquad (3.1)</math></center>

onde A, B e C são constantes, é chamada de equação quasilinear.
Um método implícito no tempo e estável numericamente de solucionar este tipo de equação foi proposto na metade do século XX por John Crank e Phyllis Nicolson.
O método de Crank-Nicolson é de segunda ordem no tempo e no espaço e é baseado em diferenças centradas no espaço e regra trapezoidal no tempo.
Como provaremos mais adiante, este método é incondicionalmente estável.

A equação diferencial finita de Nagumo é dada pela segunda derivada espacial obtida a partir de uma combinação da derivada temporal no passo <math> n</math> e <math> n+1</math>:

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} = (\theta) \frac{u_{i+1}^n -2u_i^{n} + u_{i+1}^n}{\Delta x^2} + (1-\theta) \frac{u_{i+1}^{n+1} -2u_i^{n+1} + u_{i-1}^{n+1}}{\Delta x^2} \qquad(3.2)</math></center>

A partir da equação (3.2), podemos obter os diferentes métodos de diferenciação da equação de Nagumo ajustando o valor de <math>\theta</math>.
*<math>\theta = 1</math> é o método explícito.
*<math>\theta = 0</math> é o método implícito.
*<math>\theta = \frac{1}{2}</math> é o método Crank-Nicolson.

Substituindo o valor de <math>\theta</math> na equação (3.2), temos:

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} = \frac{u_{i+1} ^n -2u_i ^{n} + u_{i+1} ^n}{\Delta x^2} + u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.3)</math></center>

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} = \frac{u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1}}{\Delta x^2} + u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha) \qquad(3.4)</math></center>

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} =\frac{u_{i+1} ^n -2u_i ^{n} + u_{i+1} ^n}{2\Delta x^2} + \frac{u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1}}{2\Delta x^2} + u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.5)</math></center>

Para nosso problema, usaremos o método pseudo Crank-Nicolson, devido ao termo não-linear presente em (3.5).

=== Análise de Estabilidade de Von Neumann ===

Ao analisar a estabilidade do método, precisamos verificar se o erro gerado no procedimento não aumenta quando testamos os extremos da solução.
Para realizar esta verificação, vamos utilizar o método de Von Neumann.
A primeira coisa a se fazer é linearizar o termo não linear da equação (3.5):

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} =\frac{u_{i+1} ^n -2u_i ^{n} + u_{i+1} ^n}{2\Delta x^2} + \frac{u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1}}{2\Delta x^2} -\alpha u_{i} ^{n}\qquad(3.6)</math></center>

Se substituírmos <math>u_{i} ^{n} = \Psi (t)e^{j\gamma \Delta x}</math></div> no tempo <math>t</math>, onde <math>j = \sqrt{(-1)}</math> e <math>\gamma > 0</math> o método se mostra estável pelo seguinte:

<center><math>\bigg| \frac{\Psi (t + \Delta t)}{\Psi (t)} \bigg| \le 1 </math></center>

Substituindo <math>u_{i} ^{n} = \Psi (t)e^{j\gamma x} </math> na equação (3.6), temos:

<center><math>\begin{matrix} \frac{\Psi (t + \Delta t) e^{j\gamma x}}{\Delta t}& =&\frac{\Psi (t) e^{j\gamma (x + \Delta x)} -2 \Psi (t) e^{j\gamma x} + \Psi (t) e^{j\gamma (x - \Delta x)}}{2\Delta x^2} \\ \ &+ &\frac{\Psi (t + \Delta t) e^{j\gamma (x + \Delta x)} -2\Psi (t + \Delta t)e^{j\gamma x} + \Psi (t + \Delta t) e^{j\gamma (x + \Delta x)}}{2\Delta x^2} -\alpha \Psi (t) e^{j\gamma x}\end{matrix}\qquad(3.7)</math></center>

Substituindo <math>r = {\Delta t}/{\Delta x^2}</math> em (3.7), temos:

<center><math>\begin{matrix} \Psi (t + \Delta t) e^{j\gamma x} &=& r \bigg(\frac{\Psi (t) e^{j\gamma (x + \Delta x)} -2\Psi (t) e^{j\gamma x} + \Psi (t) e^{j\gamma (x - \Delta x)}}{2}\bigg) \\ \ &
+& r\bigg( \frac{\Psi (t + \Delta t) e^{j\gamma (x + \Delta x)} -2\Psi (t + \Delta t)e^{j\gamma x} + \Psi (t + \Delta t) e^{j\gamma (x + \Delta x)}}{2}\bigg) -\Delta t \alpha \Psi (t) e^{j\gamma x}\end{matrix}\qquad(3.8)</math></center>

Podemos notar que há <math>e^{j\gamma x} </math> em todos os termos. Simplificando-os e movendo todos os termos com <math>\Psi (t + \Delta t) </math> para o lado esquerdo e todos com <math>\Psi (t) </math> para a direita, a equação (3.8) toma a forma:

<center><math>\Psi (t + \Delta t) - r \bigg(\frac{\Psi (t + \Delta t)e^{j\gamma \Delta x} -2\Psi (t + \Delta t) + \Psi (t + \Delta t)e^{-j\gamma \Delta x}}{2}\bigg) = r \bigg( \frac{\Psi (t) e^{j\gamma \Delta x} -2\Psi (t) + \Psi (t) e^{-j\gamma \Delta x}}{2}\bigg) + \Psi (t) -\Delta t \alpha \Psi (t)\qquad(3.9)</math></center>

Fatorando os termos, ficamos com:

<center><math>\Psi (t + \Delta t) \bigg( 1 - r \bigg( \frac{e^{j\gamma \Delta x} -2 + e^{-j\gamma \Delta x}}{2}\bigg) \bigg) = \Psi (t) \bigg( r \bigg( \frac{e^{j\gamma \Delta x} - 2 + e^{-j\gamma \Delta x}}{2} \bigg) + 1 -\Delta t \alpha \bigg) \qquad(3.10)</math></center>

Utilizando a identidade trigonométrica <math> cos(x) = \frac{e^{jx} + e^{-jx}}{2}</math>, ficamos com:

<center><math>\Psi (t + \Delta t) (1 - r (cos(\gamma \Delta x) -1)) = \Psi (t) (r (cos(j\gamma \Delta x) -1 ) +1 -\Delta t \alpha)\qquad(3.11)</math></center>

Aplicando <math>1 - cos(2u) = 2 sin^2 (u)</math>, temos:

<center><math>\Psi (t + \Delta t) \bigg(1 + 2r sin^2 \bigg( \frac{\gamma \Delta x}{2}\bigg) \bigg) = \Psi (t) \bigg(1 - 2r sin^2 \bigg( \frac{\gamma \Delta x}{2}\bigg) -\Delta t \alpha\bigg)\qquad(3.12)</math></center>

Pode-se demonstrar que, a partir da equação (3.12), tem-se:

<center><math>\frac{\Psi (t + \Delta t) e^{j\gamma x}}{\Delta t} =\frac{1 - 2r sin^2( \frac{\gamma \Delta x}{2}) -\Delta t \alpha}{1 + 2r sin^2( \frac{\gamma \Delta x}{2})}\qquad(3.13)</math></center>

Como:

<center><math>\bigg| \frac{1 - 2r sin^2( \frac{\gamma \Delta x}{2}) -\Delta t \alpha}{1 + 2r sin^2( \frac{\gamma \Delta x}{2})}\bigg| \le 1 \qquad(3.13)</math></center>

prova-se que o método Crank-Nicolson é estável para <math>\Delta t </math> suficientemente pequenos.

=== Implementação para a equação de FitzHung-Nagumo ===

Para a implementação na equação de FitzHung-Nagumo, as seguintes derivadas parciais são utilizadas:

<center><math>u(x_{i} , t_{n}) = u_{i}^n \qquad(3.14)</math></center>

<center><math>\frac{\partial u}{\partial t} (x_{i} , t_{n})= \frac{u_{i}^{n+1}- u_{i}^{n}}{\Delta t}\qquad(3.15)</math></center>

<center><math>\frac{\partial^2 u}{\partial^2 x} (x_{i} , t_{n})= \frac{1}{2}\bigg(\frac{u_{i+1}^{n+1}- 2u_{i}^{n+1}+u_{i-1}^{n+1}}{\Delta x^2} + \frac{u_{i+1}^{n}- 2u_{i}^{n}+u_{i-1}^{n}}{\Delta x^2}\bigg)\qquad(3.16)</math></center>

Neste ponto, utilizamos uma aproximação para a derivada temporal e a derivada espacial de segunda ordem e substituímos as equações (3.14), (3.15) e (3.16) na equação (2).

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} =\frac{1}{2\Delta x^2}(u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1})+ \frac{1}{2\Delta x^2}(u_{i+1} ^{n} -2u_i ^{n} + u_{i-1} ^{n})
u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.17)</math></center>

Multiplicando por <math>\Delta t</math>, temos:

<center><math>u_i ^{n+1} - u_i ^{n} =\frac{\Delta t}{2\Delta x^2}(u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1})+ \frac{\Delta t}{2\Delta x^2}(u_{i+1} ^{n} -2u_i ^{n} + u_{i-1} ^{n})
\Delta t u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.18)</math></center>

Fazendo a substituição:

<center><math> r = \frac{2\Delta t}{\Delta x^2}</math></center>

e rearranjando, temos:

<center><math>-ru_{i+1}^{n+1} + (1+2r)u_{i}^{n+1} - ru_{i-1}^{n+1} = r(u_{i+1}^{n} + u_{i-1}^{n}) + (1-2r)u_{i}^{n} + \Delta tu_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.19)</math></center>

Podemos, então representar a equação (3.19) em forma matricial:

<center><math>\mathbf{AU} = \mathbf{F}</math></center>

Onde o lado esquerdo da equação é a matriz tridimensional a seguir
<center><math>
A =
\begin{bmatrix}
1+2r & -r & 0 & \cdots & 0 \\
-r & 1+2r & -r & \cdots & 0 \\
0 & -r & 1+2r & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & -r \\
0 & 0 & 0 & -r & 1+2r \\
\end{bmatrix}
</math></center>

e o lado direito:
<center><math>
F =
\begin{bmatrix}
0 \\
\vdots \\
r(u_{n+1 ^n} + u_{i-1 ^n}) + (1-2r)u_i ^n + \Delta t u_i ^n(1-u_i ^n)(u_i ^n - \alpha) \\
\vdots \\
1\\
\end{bmatrix}
</math></center>

Com isso, podemos encontrar as soluções da equação resolvendo:

<center><math>\mathbf{U} = \mathbf{F/A}</math></center>

===FTCS Explícito===
No método explícito, a discretização é

<math>
\begin{cases}
\frac{\partial f}{\partial t} \rightarrow \frac{f(x, t+\Delta t) - f(x,t)}{\Delta t}\\
\frac{\partial^2 f}{\partial x^2} \rightarrow \frac{f(x, t+\Delta t) - 2f(x,t) + f(x-\Delta t}{(\Delta x)^2}
\end{cases}
</math>

Substituindo na EDU e utilizando a notação de índices (<math>n</math> é índice temporal e <math>j</math> é o índice temporal), temos

<math>
\frac{f_j^{n+1} - f_j^n}{\Delta t} = D\frac{f_{j+1}^n - 2f_j^n + f_{j-1}^n}{(\Delta x)^2}
</math>

Agrupando <math>\Delta t,\ \Delta x\ e\ D</math> em uma constante <math>k</math> e deixando o passo temporal futuro em função do passo anterior explicitamente, temos

<math>
f_j^{n+1} = f_j^n+k(f_{j+1}^n - 2f_j^n + f_{j-1}^n).
</math>

Utilizando esta equação é possível calcular o valor de <math>f</math> iterativamente para todo intervalo de tempo.

===FTCS Implícito===
No método implícito, também chamado de BTCS (para trás no tempo centrado no espaço, do inglês: Backwards Time Centered Space), pois o passo temporal é dado para trás, assim a derivada temporal é

<math>
\begin{cases}
\frac{\partial f}{\partial t} \rightarrow \frac{f_j^{n} - f_j^{n-1}}{\Delta t}\\
\end{cases}
</math>

Resultando na seguinte EDU

<math>
f_j^{n} = f_j^{n-1}+k(f_{j+1}^n - 2f_j^n + f_{j-1}^n).
</math>

Reescrevendo os índices <math>n</math> temos

<math>
f_j^{n+1} = f_j^{n}+k(f_{j+1}^{n+1} - 2f_j^{n+1} + f_{j-1}^{n+1})
</math>

Onde <math>f_j^{n+1}</math> está escrita implicitamente, e pode ser calculada utilizando o Algoritmo de Thomas [https://moodle.ufrgs.br/pluginfile.php/3653465/mod_resource/content/1/thomas.pdf].

===Equação de Recuperação 1D===

No modelo de Fitzhugh-Nagumo, utilizou-se o mesmo algoritmo de Crank-Nicolson utilizado na Equação de Nagumo, na seção 3.1. Entretanto, ainda é necessário resolver a equação diferencial de <math>w</math>, dada por

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

Para isso utilizou-se o método FTCS explícito. Discretizando a derivada, temos

<math>
\frac{\partial w}{\partial t} \rightarrow \frac{w(x,t + \Delta t)-w(x,t)}{\Delta t}
</math>

Utilizando a notação de índices e aplicando na equação da recuperação, temos

<math>
\frac{w_j^{n+1}-w_j^n}{\Delta t} = \epsilon(v_j^n - \gamma w_j^n)
</math>

<math>
w_j^{n+1}-w_j^n = \Delta t \epsilon(v_j^n - \gamma w_j^n).
</math>

Assim, uma vez conhecido <math>v_j^n</math>, é possível calcular <math>w</math> para qualquer intervalo de tempo.

Vamos verificar as condições de estabilidade do método utilizando o os modos de Fourier:
*substituindo <math>w_j^n=A^ne^{i(jq\Delta x)},</math> onde <math>i</math> é a unidade imaginária;
*desconsiderando <math>v_j^n</math>,
temos

<math>
A^{n+1}e^{i(jq \Delta x)} = -\Delta t \epsilon \gamma e^{i(lq \delta x)}
</math>

Dividindo ambos os lados por <math>A^ne^{i(jq\Delta x)},\ </math> temos

<math>
\frac{A^{n+1}}{A^n} = -\Delta t \epsilon \gamma
</math>

Para que o método seja estável, temos que <math>\Bigg|\frac{A^{n+1}}{A^n}\Bigg| < 1,\ </math> logo

<math>
\Bigg|\frac{A^{n+1}}{A^n}\Bigg| =\Bigg|-\Delta t \epsilon \gamma \Bigg|< 1.
</math>

Assim, temos que a condição de estabilidade sobre o método é tal que

<math>
\frac{-1}{\epsilon \gamma}<\Delta t < \frac{1}{\epsilon \gamma}.
</math>

===Modelo FitzHung-Nagumo 2D===
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 de Nagumo em 2D. Assumindo <math>\Delta x = \Delta y</math>, temos

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

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

<math>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} - 4 v_{l,j}^{n}\right] + \Delta t f_{N l,j}^{n} - \Delta t w_{l,j}^{n}</math>

onde os índices <math>l</math> e <math>j</math> se referem às coordenadas <math>x</math> e <math>y</math>, respectivamente, e o índice <math>n</math> refere-se ao tempo.

Vamos verificar as condições de estabilidade do método utilizando os modos de Fourier:
*substituindo <math>v_{l,j}^{n} = A^{n}e^{i(xq_x +yq_y)} = A^{n}e^{i\beta}</math>
*usando <math>r = \frac{D \Delta t}{(\Delta x)^2}</math>
*por simplicidade, como fizemos na seção 3.1, vamos linearizar <math>f_{N l,j}^{n} \rightarrow -\alpha v_{l,j}^{n} =-\alpha A^{n}e^{i\beta}</math>
*por fim, desconsiderando <math>\Delta t w_{l,j}^{n}</math>, temos

<math>A^{n+1}e^{i\beta} = A^{n}e^{i\beta} + rA^{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)} - 4 e^{i\beta}\right] -\Delta t \alpha A^{n}e^{i\beta}</math>

<math>A^{n+1}e^{i\beta} = A^{n}e^{i\beta} + rA^{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}</math>

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

<math>\frac{A^{n+1}}{A^{n}} = 1 + r\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</math>

Sabendo as seguintes identidades trigonométricas:

*<math>e^{-i\theta} + e^{i\theta} = 2\cos(\theta)</math>
*<math>\cos(2\theta) = 1 - 2\sin^2(\theta)</math>

e aplicando à equação, obtemos

<math>\frac{A^{n+1}}{A^{n}} = 1 + 2r\left[ \cos(\Delta xq_x) + \cos(\Delta yq_y) - 2\right] -\Delta t \alpha</math>

<math>\frac{A^{n+1}}{A^{n}} = 1 + 2r\left[ (1 - 2\sin^2(\Delta xq_x/2)) + (1 - 2\sin^2(\Delta yq_y/2)) - 2\right] -\Delta t \alpha</math>

<math>\frac{A^{n+1}}{A^{n}} = 1 + 4r\left[ \sin^2(\Delta xq_x/2) + \sin^2(\Delta yq_y/2)\right] -\Delta t \alpha.</math>

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

<math>\left|1 -\Delta t \alpha + 4r\left[ \sin^2(\Delta xq_x/2) + \sin^2(\Delta yq_y/2)\right]\right| < 1</math>.

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

<math>\frac{\alpha \Delta t}{8} < r < \frac{2 + \alpha \Delta t}{8}</math>.

==Resultados==
===Equação de Nagumo 1D===
(explicar os testes e os gráficos/animações --> Natália)

[[Arquivo:N_fino_abaixolimiar.gif|800px|thumb|center|'''Figura 4''' — estímulo <math>\ell = 20</math> <math>I = 0.05</math>]]

[[Arquivo:N_fino_acimalimiar.gif|800px|thumb|center|'''Figura 5''' — estímulo <math>\ell = 20</math> <math>I = 0.2</math>]]

[[Arquivo:N_fino_iguallimiar.gif|800px|thumb|center|'''Figura 6''' — estímulo <math>\ell = 20</math> <math>I = 0.1</math>]]

[[Arquivo:N_largo_abaixolimiar.gif|800px|thumb|center|'''Figura 7''' — estímulo <math>\ell = 100</math> <math>I = 0.05</math>]]

[[Arquivo:N_largo_acimalimiar.gif|800px|thumb|center|'''Figura 8''' — estímulo <math>\ell = 100</math> <math>I = 0.2</math>]]

[[Arquivo:N_largo_iguallimiar.gif|800px|thumb|center|'''Figura 9''' — estímulo <math>\ell = 100</math> <math>I = 0.1</math>]]

===Modelo FitzHung-Nagumo 1D===
Para a implementação da simulação do modelo FN em 1D utilizou-se <math>\Delta t = 0.1</math> (com <math>0<t<30000</math>), <math>\Delta x = 1</math> (com <math>0<x<120</math>) e os demais parâmetros foram idênticos àqueles encontrados na ref.[1] <ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </ref>. Eles são

<math>
\begin{cases}
\alpha = 0.1\\
\epsilon = 0.005\\
\gamma = 0.0001\\
D = 0.1\\
\end{cases}
</math>

Foram aplicadas condições de contorno nulas nas bordas do intervalo de <math>x</math>.

[[Arquivo:FN_vXt.png|800px|thumb|center|'''Figura 10''' — Evolução temporal do potencial elétrico em um ponto da membrana do axônio quando excitado por um estímulo acima do limiar de potencial.]]

Na figura 10 vemos a evolução temporal do potencial elétrico em um ponto da membrana do axônio (<math>x=15</math>) quando excitado acima do limiar de potencial. Note que a curva do PA é muito semelhante àquela da na Figura 2, apresentando vários dos aspectos característicos, como a assimetria entre os períodos de polarização (ascensão rápida) e repolarização (decaimento mais lento) e um período refratário, onde o potencial lentamente se recupera em direção do PR. Podemos assim entender a escolha dos valores dos parâmetros <math>\epsilon</math> e <math>\gamma</math>, que controlam a velocidade de crescimento da variável de recuperação (assim como se pode ver na equação (3)). Experimentalmente, sabemos que o período refratário é muito mais longo do que o período de despolarização
<ref name=QUILLFELDT> https://www.ufrgs.br/mnemoforos/arquivos/potenciais2005.pdf Jorge A. Quillfeldt,"ORIGEM DOS POTENCIAIS ELÉTRICOS DAS CÉLULAS NERVOSAS" </ref>
, logo precisamos que a intensidade da variável <math>w</math> aumente muito lentamente. Para isso, escolhe-se valores pequenos de <math>\epsilon</math> e <math>\gamma</math>, pois caso contrário o potencial da membrana celular não teria tempo de atingir o pico do PA.

A velocidade de crescimento da variável de recuperação também explica a ausência do "tudo ou nada" no PA da figura 10. Na figura 11, vê-se a evolução temporal de um PA com <math>\epsilon = 0.001</math> (5 vezes menor do que aquele da figura 10).

[[Arquivo:FN_vXt_Epequeno.png|800px|thumb|center|'''Figura 11''' — Evolução temporal do potencial elétrico em um ponto da membrana do axônio para <math>\epsilon = 0.001</math>.]]

Na figura 11, como a variável de recuperação cresce muito mais lentamente, além de demorar mais tempo para a membrana celular sair do período refratário, vemos que o pico do PA de ação é muito mais próximo do valor máximo. Isso ocorre porque a variável de recuperação cresce mais lentamente, possibilitando que a variável rápida atinja um valor maior antes de começar a diminuir devido à ação de <math>w</math>.

[[Arquivo:FN_abaixolim_apos.gif|800px|thumb|center|'''Figura 12''' — Simulação de um estímulo excitando um axônio, segundo o modelo de Fitzhugh-Nagumo. Essa simulação foi feita para <math>\alpha = 0.1</math> e um estímulo inicial de largura <math>\ell = 10</math> e intensidade <math>I=0.05</math> (abaixo do limiar). Note que o estímulo não desencadeia um PA e logo desaparece.]]

[[Arquivo:FN_acimalim_apos.gif|800px|thumb|center|'''Figura 13''' — Simulação de um estímulo excitando um axônio, segundo o modelo de Fitzhugh-Nagumo. Essa simulação foi feita para <math>\alpha = 0.1</math> e um estímulo inicial de largura <math>\ell = 10</math> e intensidade <math>I=0.2</math> (acima do limiar). Note que quando o estímulo está acima do limiar, ele gera um PA que se propaga ao longo do axônio.]]

Nas animações das figuras 12 e 13 vemos a evolução espacial do potencial elétrico com o decorrer do tempo.
Na figura 12, o estímulo inicial está abaixo do limiar de voltagem, desta forma ele não gera um PA e a célula volta rapidamente para o seu PR. Na figura 13 o estímulo inicial está acima do limiar, logo um PA é gerado e o pulso elétrico é propagado ao longo do axônio.

===Modelo FitzHung-Nagumo 2D===

As simulação 2D tinham por objetivo simular a propagação de um estimulo inicial <math>I = 0.2</math> de <math>\ell = 5\times 50</math> de espessura. Para se aproximar da simetria do problema(axônio como um cilindro por onde o PA se propaga) aplicamos condições periódicas na vertical, e mantivemos as bordas nulas na horizontal.

Os parâmetros adotados para <math>\alpha, \epsilon, \gamma, D, \Delta x </math> e <math> \Delta t</math> foram os mesmos da secção 5.2.

[[Arquivo:PA_Sob_Membrana2D.png|800px|thumb|center|'''Figura 14''' — Imagens de simulação de PA propagando sobre membrana 2D segundo modelo de FitzHung-Nagumo]]

Como podemos ver o PA se propaga assumindo um front de alta voltagem(despolarização e repolarização) estreito, e logo após uma região mais esparsa(de roxo mais escuro) indicando a repolarização da membrana da célula.

==Discussão==

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 <ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>

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

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 <math>\alpha > 0</math> <ref name=QUILLFELDT> https://www.ufrgs.br/mnemoforos/arquivos/potenciais2005.pdf Jorge A. Quillfeldt,"ORIGEM DOS POTENCIAIS ELÉTRICOS DAS CÉLULAS NERVOSAS" </ref>,que como mencionado na seçã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 <math>\alpha < 0</math>, mas embora essa abordagem leve a uma gama maior de resultados para estímulo contínuo
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>
<ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and Biology Society, Aug 2014, Chicago, United States. ffhal-00998828</ref> o sentido fisiológico é prejudicado, principalmente pelo fato de não possuir um limiar
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>.

==Programas==

===Cranck Nicolson e FTCS Modelo FitzHung-Nagumo 1D===
<source lang="python">

import numpy as np
import matplotlib
from matplotlib import pyplot as plt
from matplotlib.animation import ArtistAnimation

#Definindo Polinômio de Nagumo
def Pol(u):
a = 0.1
return u*(1 - u)*(u - a)

#Declarando Constantes
L = 100
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)

#Definindo função que calcula a equação de recuperação 1D
def FitzHung_FTCS(v, w):
D = 0.1
dx = 1
r = (D*dt)/(dx**2)
L = len(v)
aux = w[:]
for i in range(1, L-1):
aux[i] = w[i] + (r/(D*2))*0.005*(v[i] - 0.0001*w[i])

for j in range(1, L-1):
w[j] = aux[j]

return (w)

#Criando matriz tridiagonal para o método de Thomas
for i in range(L):
for j in range(L):
if i == j:
A[i,j] = 1 - 2*r
if j == i+1:
A[i,j] = r
if j == i-1:
A[i,j] = r

#Definindo função que calcula Equação de Nagumo 1D por método de Thomas
def Nagumo_CN(d,w,A):
L = len(d)
d = np.array(d)
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)

aux = np.dot(A, d)
aux[0] = d[0]
aux[L-1] = d[L-1]
for i in range(L):
d[i] = aux[i] + (r*Pol(aux[i]))/(D*2) - (r*w[i])/(2*D) + (r/(2*D))*I[i]

a = (-1*r)*(np.ones(L, dtype = np.float32))

b = (2*r + 1)*(np.ones(L, dtype = np.float32))

c = (-1*r)*(np.ones(L, dtype = np.float32))
c[L-1] = 0
g = [x for x in d]

c[0] = c[0]/b[0]
for i in range(1,L-1) :
c[i] = c[i]/(b[i] - (c[i-1]*a[i]))

for i in range(1,L-1) :
d[i] = (d[i] - (d[i-1]*a[i]))/(b[i] - (c[i-1]*a[i]))

g[L-1] = d[L-1]
for i in range(L-2,0,-1) :
g[i] = (d[i] - c[i]*g[i+1])

return (g)

#calculando o modelo FN
f = np.zeros(L, dtype = np.float32)
w = np.zeros(L, dtype = np.float32)
f[1:11] = 0.2

h = np.zeros(30000, dtype = np.float32)

T = np.arange(30000, dtype = np.float32)

x = np.arange(L, dtype = np.float32)
t = 0
h[0] = f[15]
for m in range(1, 30000):

aux = f[:]
f = Nagumo_CN(f, w, A)
w = FitzHung_FTCS(aux , w)
h[m] = f[15]
t = m*dt

#plot das figuras 10 e 11

ones = np.ones(30000)
zeros = np.zeros(30000)
plt.title('Modelo FN, $\\upsilon \\times t$ ($\ell = 10$, $I = 0.2$)')
plt.xlabel('t')
plt.ylabel('$\\upsilon$')
plt.plot(T, h)
plt.plot(T, ones,'r--', linewidth=0.8)
plt.plot(T, zeros, 'k--', linewidth = 0.8)
plt.savefig('nomedoarquivo.png')

#plot das animações

f = np.zeros(L, dtype = np.float32)
w = np.zeros(L, dtype = np.float32)
f[1:11] = 0.2
x = np.arange(L, dtype = np.float32)
t = 0
n = 1

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
frames = []
plt.title('Model FN, $\\ell$ = 10, I = 0.2')
plt.ylabel('$\\upsilon$')
plt.xlabel('x')
for m in range(1, 10000):

if (t == n):
curve = ax.plot(x, f, 'c')
frames.append(curve)
n = n + 1
aux = f[:]
f = Nagumo_CN(f, w, A)
w = FitzHung_FTCS(aux , w)

t = m*dt

animacao = ArtistAnimation(fig, frames, interval=50, blit=True)
animacao.save('nomedoarquivo.mp4')

</source>

===FTCS Modelo FitzHung-Nagumo 2D===
<source lang="python">

import numpy as np
import matplotlib
from matplotlib import pyplot as plt
import seaborn as sb

#Definindo Polinômio de Nagumo
def Pol(u):
a = 0.1
return u*(1 - u)*(u - a)

#Declarando Constantes
L = 100
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)

#Definindo função que calcula a equação de recuperação 2D
def FitzHung2D_FTCS(v, w):
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)
L = len(v)

aux = w[:]
for i in range(1, L-1):
for j in range(1, L-1):
aux[i,j] = w[i,j] + dt*0.005*(v[i,j] - 0.0001*w[i,j])

for i in range(1, L-1):
for j in range(1, L-1):
w[i,j] = aux[i,j]

return (w)

#Definindo função que calcula a equação de Nagumo em 2D como condições periódicas em y
def Nagumo2D_FTCS(v, w):
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)
L = len(v)

aux = v[:,:]
for j in range(1,L-1):
aux[0,j] = v[0,j] + r*(v[ L-1 ,j] + v[1,j] + v[0,j-1] + v[0,j+1] -4*v[0,j]) + dt*Pol(v[0,j]) - dt*w[0,j]

aux[0,L-1] = v[0,L-1] + r*(v[ L-1 ,L-1] + v[1,L-1] + v[0,L-2] + 0 -4*v[0,L-1]) + dt*Pol(v[0,L-1]) - dt*w[0,L-1]

for i in range(1, L-1):
for j in range(1, L-1):
aux[i,j] = v[i,j] + r*(v[i-1,j] + v[i+1,j] + v[i,j-1] + v[i,j+1] -4*v[i,j]) + dt*Pol(v[i,j]) - dt*w[i,j]

for j in range(1,L-1):
aux[L-1,j] = v[L-1,j] + r*(v[L-2,j] + v[0,j] + v[L-1,j-1] + v[L-1,j+1] -4*v[L-1,j]) + dt*Pol(v[L-1,j]) - dt*w[L-1,j]

aux[L-1,L-1] = v[L-1,L-1] + r*(v[L-2,L-1] + v[0,L-1] + v[L-1,L-2] + 0 -4*v[L-1,L-1]) + dt*Pol(v[L-1,L-1]) - dt*w[L-1,L-1]

for i in range(1, L):
for j in range(1, L-1):
v[i,j] = aux[i,j]

return (v)

#Função criada para evitar o uso de estruturas booleanas dentro do loop, que torna a simulação mais demorada
def Para_Evitar_Ifs(g):
intervalo = 20
T = g[0]
v = g[1]
w = g[2]
for t in range(T,T + intervalo):
aux = v[:,:]
v = Nagumo2D_FTCS(v, w)
w = FitzHung2D_FTCS(aux , w)

return ([T+intervalo,v,w])

#Criando os espaço e o grid de fotos
x = np.arange(L, dtype = np.float32)
y = np.arange(L, dtype = np.float32)
T = 0
dl = 4
p1 = (dl/2) - 1
p2 = (dl/2) + 1
tamanho = 2*L/dl

x,y = np.meshgrid(x,y)

f, (ax1,ax2,ax3,ax4) = plt.subplots(1,4,sharey=True, figsize=(24,4))
f.suptitle('Pulso Lateral de Tamanho '+str(int(tamanho))+'')

#Declarando condições iniciais e condições de contorno não periódicas
v = np.zeros((L,L))
v[int(p1*tamanho/2)-6:int(p2*tamanho/2),1:5] = 0.2

w = np.zeros((L,L))

#Plotando os 4 primeiros mapas de calor
g = [T,v,w]
graph1 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax1)
graph1.set_title('t='+str(g[0])+'')
graph1.invert_yaxis()

g = Para_Evitar_Ifs(g)
graph2 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax2)
graph2.set_title('t='+str(g[0])+'')
graph2.invert_yaxis()

g = Para_Evitar_Ifs(g)
graph3 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax3)
graph3.set_title('t='+str(g[0])+'')
graph3.invert_yaxis()

g = Para_Evitar_Ifs(g)
graph4 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax4)
graph4.set_title('t='+str(g[0])+'')
graph4.invert_yaxis()

plt.savefig('Lateral_1.png')
plt.show()

</source>

==Referências==
<references/>

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

2021-04-04T04:00:45Z

Nataliaferrazzo:

'''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)<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| A simple spiking neuron model: sodium channels alone </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 Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </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>
<ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </ref>
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>

<center><math>
\begin{cases}
\frac{\partial v}{\partial t} = D\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)
</math></center>

Onde <math>0<\alpha<1</math> é o limiar de potencial, <math>D</math> é uma constante de difusão, <math>\epsilon</math> e <math>\gamma</math> são parâmetros positivos a serem ajustados —relacionados à velocidade de atuação da variável de recuperação— e <math>I</math> é a magnitude do estímulo externo <ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </ref>
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>

== Método Crank-Nicolson ==

Uma equação diferencial parcial da forma
<center><math>Au_{xx} + Bu_{xy} + Cu_{yy} = f(x, y, u_x, u_y) \qquad (3.1)</math></center>

onde A, B e C são constantes, é chamada de equação quasilinear.
Um método implícito no tempo e estável numericamente de solucionar este tipo de equação foi proposto na metade do século XX por John Crank e Phyllis Nicolson[https://pt.wikipedia.org/wiki/M%C3%A9todo_de_Crank%E2%80%93Nicolson#cite_note-2].
O método de Crank-Nicolson é de segunda ordem no tempo e no espaço e é baseado em diferenças centradas no espaço e regra trapezoidal no tempo[https://pt.wikipedia.org/w/index.php?title=Regra_Trapezoidal&action=edit&redlink=].
Como provaremos mais adiante, este método é incondicionalmente estável.

A equação diferencial finita de Nagumo é dada pela segunda derivada espacial obtida a partir de uma combinação da derivada temporal no passo <math> n</math> e <math> n+1</math>:

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} = (\theta) \frac{u_{i+1}^n -2u_i^{n} + u_{i+1}^n}{(\Delta x^2)} + (1-\theta) \frac{u_{i+1}^{n+1} -2u_i^{n+1} + u_{i-1}^{n+1}}{(\Delta x^2)} \qquad(3.2)</math></center>

A partir da equação (3.2), podemos obter os diferentes métodos de diferenciação da equação de Nagumo ajustando o valor de <math>\theta</math>.
*<math>\theta = 1</math> é o método explícito.
*<math>\theta = 0</math> é o método implícito.
*<math>\theta = \frac{1}{2}</math> é o método Crank-Nicolson.

Substituindo o valor de <math>\theta</math> na equação (3.2), temos:

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} = \frac{u_{i+1} ^n -2u_i ^{n} + u_{i+1} ^n}{(\Delta x^2)} + u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.3)</math></center>

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} = \frac{u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1}}{(\Delta x^2)} + u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha) \qquad(3.4)</math></center>

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} =\frac{u_{i+1} ^n -2u_i ^{n} + u_{i+1} ^n}{(2\Delta x^2)} + \frac{u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1}}{(2\Delta x^2)} + u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.5)</math></center>

Para nosso problema, usaremos o método pseudo Crank-Nicolson, devido ao termo não-linear presente em (3.5).

=== Implementação para a equação de FitzHung-Nagumo ===

Para a implementação na equação de FitzHung-Nagumo, as seguintes derivadas parciais são utilizadas:

<center><math>u(x_{i} , t_{n}) = u_{i}^n \qquad(3.6)</math></center>

<center><math>\frac{\partial u}{\partial t} (x_{i} , t_{n})= \frac{u_{i}^{n+1}- u_{i}^{n}}{\Delta t}\qquad(3.7)</math></center>

<center><math>\frac{\partial^2 u}{\partial^2 x} (x_{i} , t_{n})= \frac{1}{2}(\frac{u_{i+1}^{n+1}- 2u_{i}^{n+1}+u_{i-1}^{n+1}}{\Delta x^2} + \frac{u_{i+1}^{n}- 2u_{i}^{n}+u_{i-1}^{n}}{\Delta x^2})\qquad(3.8)</math></center>

Neste ponto, utilizamos uma aproximação para a derivada temporal e a derivada espacial de segunda ordem e substituímos as equações (3.6), (3.7) e (2.8) na equação (2).

<center><math>\frac{u_i ^{n+1} - u_i ^{n}}{\Delta t} =\frac{1}{(2\Delta x^2)}(u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1})+ \frac{1}{(2\Delta x^2)}(u_{i+1} ^{n} -2u_i ^{n} + u_{i-1} ^{n})
u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.9)</math></center>

Multiplicando por <math>\Delta t</math>, temos:

<center><math>u_i ^{n+1} - u_i ^{n} =\frac{\Delta t}{(2\Delta x^2)}(u_{i+1} ^{n+1} -2u_i ^{n+1} + u_{i-1} ^{n+1})+ \frac{\Delta t}{(2\Delta x^2)}(u_{i+1} ^{n} -2u_i ^{n} + u_{i-1} ^{n})
\Delta t u_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.10)</math></center>

Fazendo a substituição:

<center><math> r = \frac{2\Delta t}{\Delta x^2}</math></center>

e rearranjando, temos:

<center><math>-ru_{i+1}^{n+1} + (1+2r)u_{i}^{n+1} - ru_{i-1}^{n+1} = r(u_{i+1}^{n} + u_{i-1}^{n}) + (1-2r)u_{i}^{n} + \Delta tu_{i} ^{n} (1-u_{i} ^{n})(u_{i} ^{n} - \alpha)\qquad(3.11)</math></center>

Podemos, então representar a equação (3.11) em forma matricial:

<center><math>\mathbf{AU} = \mathbf{F}</math></center>

Onde o lado esquerdo da equação é a matriz tridimensional a seguir
<center><math>
A =
\begin{bmatrix}
1+2r & -r & 0 & \cdots & 0 \\
-r & 1+2r & -r & \cdots & 0 \\
0 & -r & 1+2r & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & -r \\
0 & 0 & 0 & -r & 1+2r \\
\end{bmatrix}
</math></center>

e o lado direito:
<center><math>
F =
\begin{bmatrix}
0 \\
\vdots \\
r(u_{n+1 ^n} + u_{i-1 ^n}) + (1-2r)u_i ^n + \Delta t u_i ^n(1-u_i ^n)(u_i ^n - \alpha) \\
\vdots \\
1\\
\end{bmatrix}
</math></center>

Com isso, podemos encontrar as soluções da equação resolvendo:

<center><math>\mathbf{U} = \mathbf{F/A}</math></center>

==Método FTCS==
O FTCS (Adiante no tempo, centrado no espaço, do inglês: Forward Time Centered Space) é um método numérico para se resolver equações diferenciais parciais. Ele se baseia na discretização das derivadas parciais da equação e pode ser tanto um método explícito quanto um método implícito. Como exemplo, vamos utilizar a equação abaixo, conhecida como a equação da difusão unidimensional (EDU)

<math>\frac{\partial f}{\partial t}=D\frac{\partial^2 f}{\partial x^2}</math>

===FTCS Explícito===
No método explícito, a discretização é

<math>
\begin{cases}
\frac{\partial f}{\partial t} \rightarrow \frac{f(x, t+\Delta t) - f(x,t)}{\Delta t}\\
\frac{\partial^2 f}{\partial x^2} \rightarrow \frac{f(x, t+\Delta t) - 2f(x,t) + f(x-\Delta t}{(\Delta x)^2}
\end{cases}
</math>

Substituindo na EDU e utilizando a notação de índices (<math>n</math> é índice temporal e <math>j</math> é o índice temporal), temos

<math>
\frac{f_j^{n+1} - f_j^n}{\Delta t} = D\frac{f_{j+1}^n - 2f_j^n + f_{j-1}^n}{(\Delta x)^2}
</math>

Agrupando <math>\Delta t,\ \Delta x\ e\ D</math> em uma constante <math>k</math> e deixando o passo temporal futuro em função do passo anterior explicitamente, temos

<math>
f_j^{n+1} = f_j^n+k(f_{j+1}^n - 2f_j^n + f_{j-1}^n).
</math>

Utilizando esta equação é possível calcular o valor de <math>f</math> iterativamente para todo intervalo de tempo.

===FTCS Implícito===
No método implícito, também chamado de BTCS (para trás no tempo centrado no espaço, do inglês: Backwards Time Centered Space), pois o passo temporal é dado para trás, assim a derivada temporal é

<math>
\begin{cases}
\frac{\partial f}{\partial t} \rightarrow \frac{f_j^{n} - f_j^{n-1}}{\Delta t}\\
\end{cases}
</math>

Resultando na seguinte EDU

<math>
f_j^{n} = f_j^{n-1}+k(f_{j+1}^n - 2f_j^n + f_{j-1}^n).
</math>

Reescrevendo os índices <math>n</math> temos

<math>
f_j^{n+1} = f_j^{n}+k(f_{j+1}^{n+1} - 2f_j^{n+1} + f_{j-1}^{n+1})
</math>

Onde <math>f_j^{n+1}</math> está escrita implicitamente, e pode ser calculada utilizando o Algoritmo de Thomas [https://moodle.ufrgs.br/pluginfile.php/3653465/mod_resource/content/1/thomas.pdf].

===Equação de Recuperação 1D===

No modelo de Fitzhugh-Nagumo, utilizou-se o mesmo algoritmo de Crank-Nicolson utilizado na Equação de Nagumo, na seção 3.1. Entretanto, ainda é necessário resolver a equação diferencial de <math>w</math>, dada por

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

Para isso utilizou-se o método FTCS explícito. Discretizando a derivada, temos

<math>
\frac{\partial w}{\partial t} \rightarrow \frac{w(x,t + \Delta t)-w(x,t)}{\Delta t}
</math>

Utilizando a notação de índices e aplicando na equação da recuperação, temos

<math>
\frac{w_j^{n+1}-w_j^n}{\Delta t} = \epsilon(v_j^n - \gamma w_j^n)
</math>

<math>
w_j^{n+1}-w_j^n = \Delta t \epsilon(v_j^n - \gamma w_j^n).
</math>

Assim, uma vez conhecido <math>v_j^n</math>, é possível calcular <math>w</math> para qualquer intervalo de tempo.

Vamos verificar as condições de estabilidade do método utilizando o os modos de Fourier:
*substituindo <math>w_j^n=A^ne^{i(jq\Delta x)},</math> onde <math>i</math> é a unidade imaginária;
*desconsiderando <math>v_j^n</math>,
temos

<math>
A^{n+1}e^{i(jq \Delta x)} = -\Delta t \epsilon \gamma e^{i(lq \delta x)}
</math>

Dividindo ambos os lados por <math>A^ne^{i(jq\Delta x)},\ </math> temos

<math>
\frac{A^{n+1}}{A^n} = -\Delta t \epsilon \gamma
</math>

Para que o método seja estável, temos que <math>\Bigg|\frac{A^{n+1}}{A^n}\Bigg| < 1,\ </math> logo

<math>
\Bigg|\frac{A^{n+1}}{A^n}\Bigg| =\Bigg|-\Delta t \epsilon \gamma \Bigg|< 1.
</math>

Assim, temos que a condição de estabilidade sobre o método é tal que

<math>
\frac{-1}{\epsilon \gamma}<\Delta t < \frac{1}{\epsilon \gamma}.
</math>

===Modelo FitzHung-Nagumo 2D===
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 de Nagumo em 2D. Assumindo <math>\Delta x = \Delta y</math>, temos

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

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

<math>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} - 4 v_{l,j}^{n}\right] + \Delta t f_{N l,j}^{n} - \Delta t w_{l,j}^{n}</math>

onde os índices <math>l</math> e <math>j</math> se referem às coordenadas <math>x</math> e <math>y</math>, respectivamente, e o índice <math>n</math> refere-se ao tempo.

Vamos verificar as condições de estabilidade do método utilizando os modos de Fourier:
*substituindo <math>v_{l,j}^{n} = A^{n}e^{i(xq_x +yq_y)} = A^{n}e^{i\beta}</math>
*usando <math>r = \frac{D \Delta t}{(\Delta x)^2}</math>
*por simplicidade, como fizemos na seção 3.1, vamos linearizar <math>f_{N l,j}^{n} \rightarrow -\alpha v_{l,j}^{n} =-\alpha A^{n}e^{i\beta}</math>
*por fim, desconsiderando <math>\Delta t w_{l,j}^{n}</math>, temos

<math>A^{n+1}e^{i\beta} = A^{n}e^{i\beta} + rA^{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)} - 4 e^{i\beta}\right] -\Delta t \alpha A^{n}e^{i\beta}</math>

<math>A^{n+1}e^{i\beta} = A^{n}e^{i\beta} + rA^{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}</math>

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

<math>\frac{A^{n+1}}{A^{n}} = 1 + r\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</math>

Sabendo as seguintes identidades trigonométricas:

*<math>e^{-i\theta} + e^{i\theta} = 2\cos(\theta)</math>
*<math>\cos(2\theta) = 1 - 2\sin^2(\theta)</math>

e aplicando à equação, obtemos

<math>\frac{A^{n+1}}{A^{n}} = 1 + 2r\left[ \cos(\Delta xq_x) + \cos(\Delta yq_y) - 2\right] -\Delta t \alpha</math>

<math>\frac{A^{n+1}}{A^{n}} = 1 + 2r\left[ (1 - 2\sin^2(\Delta xq_x/2)) + (1 - 2\sin^2(\Delta yq_y/2)) - 2\right] -\Delta t \alpha</math>

<math>\frac{A^{n+1}}{A^{n}} = 1 + 4r\left[ \sin^2(\Delta xq_x/2) + \sin^2(\Delta yq_y/2)\right] -\Delta t \alpha.</math>

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

<math>\left|1 -\Delta t \alpha + 4r\left[ \sin^2(\Delta xq_x/2) + \sin^2(\Delta yq_y/2)\right]\right| < 1</math>.

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

<math>\frac{\alpha \Delta t}{8} < r < \frac{2 + \alpha \Delta t}{8}</math>.

==Resultados==
===Equação de Nagumo 1D===
(explicar os testes e os gráficos/animações --> Natália)

[[Arquivo:N_fino_abaixolimiar.gif|800px|thumb|center|'''Figura 4''' — estímulo <math>\ell = 20</math> <math>I = 0.05</math>]]

[[Arquivo:N_fino_acimalimiar.gif|800px|thumb|center|'''Figura 5''' — estímulo <math>\ell = 20</math> <math>I = 0.2</math>]]

[[Arquivo:N_fino_iguallimiar.gif|800px|thumb|center|'''Figura 6''' — estímulo <math>\ell = 20</math> <math>I = 0.1</math>]]

[[Arquivo:N_largo_abaixolimiar.gif|800px|thumb|center|'''Figura 7''' — estímulo <math>\ell = 100</math> <math>I = 0.05</math>]]

[[Arquivo:N_largo_acimalimiar.gif|800px|thumb|center|'''Figura 8''' — estímulo <math>\ell = 100</math> <math>I = 0.2</math>]]

[[Arquivo:N_largo_iguallimiar.gif|800px|thumb|center|'''Figura 9''' — estímulo <math>\ell = 100</math> <math>I = 0.1</math>]]

===Modelo FitzHung-Nagumo 1D===
Para a implementação da simulação do modelo FN em 1D utilizou-se <math>\Delta t = 0.1</math> (com <math>0<t<30000</math>), <math>\Delta x = 1</math> (com <math>0<x<120</math>) e os demais parâmetros foram idênticos àqueles encontrados na ref.[1] <ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of
FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the
IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and
Biology Society, Aug 2014, Chicago, United States. ffhal-00998828 </ref>. Eles são

<math>
\begin{cases}
\alpha = 0.1\\
\epsilon = 0.005\\
\gamma = 0.0001\\
D = 0.1\\
\end{cases}
</math>

Foram aplicadas condições de contorno nulas nas bordas do intervalo de <math>x</math>.

[[Arquivo:FN_vXt.png|800px|thumb|center|'''Figura 10''' — Evolução temporal do potencial elétrico em um ponto da membrana do axônio quando excitado por um estímulo acima do limiar de potencial.]]

Na figura 10 vemos a evolução temporal do potencial elétrico em um ponto da membrana do axônio (<math>x=15</math>) quando excitado acima do limiar de potencial. Note que a curva do PA é muito semelhante àquela da na Figura 2, apresentando vários dos aspectos característicos, como a assimetria entre os períodos de polarização (ascensão rápida) e repolarização (decaimento mais lento) e um período refratário, onde o potencial lentamente se recupera em direção do PR. Podemos assim entender a escolha dos valores dos parâmetros <math>\epsilon</math> e <math>\gamma</math>, que controlam a velocidade de crescimento da variável de recuperação (assim como se pode ver na equação (3)). Experimentalmente, sabemos que o período refratário é muito mais longo do que o período de despolarização
<ref name=QUILLFELDT> https://www.ufrgs.br/mnemoforos/arquivos/potenciais2005.pdf Jorge A. Quillfeldt,"ORIGEM DOS POTENCIAIS ELÉTRICOS DAS CÉLULAS NERVOSAS" </ref>
, logo precisamos que a intensidade da variável <math>w</math> aumente muito lentamente. Para isso, escolhe-se valores pequenos de <math>\epsilon</math> e <math>\gamma</math>, pois caso contrário o potencial da membrana celular não teria tempo de atingir o pico do PA.

A velocidade de crescimento da variável de recuperação também explica a ausência do "tudo ou nada" no PA da figura 10. Na figura 11, vê-se a evolução temporal de um PA com <math>\epsilon = 0.001</math> (5 vezes menor do que aquele da figura 10).

[[Arquivo:FN_vXt_Epequeno.png|800px|thumb|center|'''Figura 11''' — Evolução temporal do potencial elétrico em um ponto da membrana do axônio para <math>\epsilon = 0.001</math>.]]

Na figura 11, como a variável de recuperação cresce muito mais lentamente, além de demorar mais tempo para a membrana celular sair do período refratário, vemos que o pico do PA de ação é muito mais próximo do valor máximo. Isso ocorre porque a variável de recuperação cresce mais lentamente, possibilitando que a variável rápida atinja um valor maior antes de começar a diminuir devido à ação de <math>w</math>.

[[Arquivo:FN_abaixolim_apos.gif|800px|thumb|center|'''Figura 12''' — Simulação de um estímulo excitando um axônio, segundo o modelo de Fitzhugh-Nagumo. Essa simulação foi feita para <math>\alpha = 0.1</math> e um estímulo inicial de largura <math>\ell = 10</math> e intensidade <math>I=0.05</math> (abaixo do limiar). Note que o estímulo não desencadeia um PA e logo desaparece.]]

[[Arquivo:FN_acimalim_apos.gif|800px|thumb|center|'''Figura 13''' — Simulação de um estímulo excitando um axônio, segundo o modelo de Fitzhugh-Nagumo. Essa simulação foi feita para <math>\alpha = 0.1</math> e um estímulo inicial de largura <math>\ell = 10</math> e intensidade <math>I=0.2</math> (acima do limiar). Note que quando o estímulo está acima do limiar, ele gera um PA que se propaga ao longo do axônio.]]

Nas animações das figuras 12 e 13 vemos a evolução espacial do potencial elétrico com o decorrer do tempo.
Na figura 12, o estímulo inicial está abaixo do limiar de voltagem, desta forma ele não gera um PA e a célula volta rapidamente para o seu PR. Na figura 13 o estímulo inicial está acima do limiar, logo um PA é gerado e o pulso elétrico é propagado ao longo do axônio.

===Modelo FitzHung-Nagumo 2D===

As simulação 2D tinham por objetivo simular a propagação de um estimulo inicial <math>I = 0.2</math> de <math>\ell = 5\times 50</math> de espessura. Para se aproximar da simetria do problema(axônio como um cilindro por onde o PA se propaga) aplicamos condições periódicas na vertical, e mantivemos as bordas nulas na horizontal.

Os parâmetros adotados para <math>\alpha, \epsilon, \gamma, D, \Delta x </math> e <math> \Delta t</math> foram os mesmos da secção 5.2.

[[Arquivo:PA_Sob_Membrana2D.png|800px|thumb|center|'''Figura 14''' — Imagens de simulação de PA propagando sobre membrana 2D segundo modelo de FitzHung-Nagumo]]

Como podemos ver o PA se propaga assumindo um front de alta voltagem(despolarização e repolarização) estreito, e logo após uma região mais esparsa(de roxo mais escuro) indicando a repolarização da membrana da célula.

==Discussão==

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 <ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>

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

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 <math>\alpha > 0</math> <ref name=QUILLFELDT> https://www.ufrgs.br/mnemoforos/arquivos/potenciais2005.pdf Jorge A. Quillfeldt,"ORIGEM DOS POTENCIAIS ELÉTRICOS DAS CÉLULAS NERVOSAS" </ref>,que como mencionado na seçã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 <math>\alpha < 0</math>, mas embora essa abordagem leve a uma gama maior de resultados para estímulo contínuo
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>
<ref name=PARAMETERS> https://hal.inria.fr/hal-00998828/document Binbin Xu, Stéphane Binczak, Sabir Jacquir, Oriol Pont, Hussein Yahia. Parameters Analysis of FitzHugh-Nagumo Model for a Reliable Simulation. 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC’14), IEEE Engineering in Medicine and Biology Society, Aug 2014, Chicago, United States. ffhal-00998828</ref> o sentido fisiológico é prejudicado, principalmente pelo fato de não possuir um limiar
<ref name=SCHOLARPEDIA> http://www.scholarpedia.org/article/FitzHugh-Nagumo_model Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349., "FITZHUGH-NAGUMO MODEL" </ref>.

==Programas==

===Cranck Nicolson e FTCS Modelo FitzHung-Nagumo 1D===
<source lang="python">

import numpy as np
import matplotlib
from matplotlib import pyplot as plt
from matplotlib.animation import ArtistAnimation

#Definindo Polinômio de Nagumo
def Pol(u):
a = 0.1
return u*(1 - u)*(u - a)

#Declarando Constantes
L = 100
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)

#Definindo função que calcula a equação de recuperação 1D
def FitzHung_FTCS(v, w):
D = 0.1
dx = 1
r = (D*dt)/(dx**2)
L = len(v)
aux = w[:]
for i in range(1, L-1):
aux[i] = w[i] + (r/(D*2))*0.005*(v[i] - 0.0001*w[i])

for j in range(1, L-1):
w[j] = aux[j]

return (w)

#Criando matriz tridiagonal para o método de Thomas
for i in range(L):
for j in range(L):
if i == j:
A[i,j] = 1 - 2*r
if j == i+1:
A[i,j] = r
if j == i-1:
A[i,j] = r

#Definindo função que calcula Equação de Nagumo 1D por método de Thomas
def Nagumo_CN(d,w,A):
L = len(d)
d = np.array(d)
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)

aux = np.dot(A, d)
aux[0] = d[0]
aux[L-1] = d[L-1]
for i in range(L):
d[i] = aux[i] + (r*Pol(aux[i]))/(D*2) - (r*w[i])/(2*D) + (r/(2*D))*I[i]

a = (-1*r)*(np.ones(L, dtype = np.float32))

b = (2*r + 1)*(np.ones(L, dtype = np.float32))

c = (-1*r)*(np.ones(L, dtype = np.float32))
c[L-1] = 0
g = [x for x in d]

c[0] = c[0]/b[0]
for i in range(1,L-1) :
c[i] = c[i]/(b[i] - (c[i-1]*a[i]))

for i in range(1,L-1) :
d[i] = (d[i] - (d[i-1]*a[i]))/(b[i] - (c[i-1]*a[i]))

g[L-1] = d[L-1]
for i in range(L-2,0,-1) :
g[i] = (d[i] - c[i]*g[i+1])

return (g)

#calculando o modelo FN
f = np.zeros(L, dtype = np.float32)
w = np.zeros(L, dtype = np.float32)
f[1:11] = 0.2

h = np.zeros(30000, dtype = np.float32)

T = np.arange(30000, dtype = np.float32)

x = np.arange(L, dtype = np.float32)
t = 0
h[0] = f[15]
for m in range(1, 30000):

aux = f[:]
f = Nagumo_CN(f, w, A)
w = FitzHung_FTCS(aux , w)
h[m] = f[15]
t = m*dt

#plot das figuras 10 e 11

ones = np.ones(30000)
zeros = np.zeros(30000)
plt.title('Modelo FN, $\\upsilon \\times t$ ($\ell = 10$, $I = 0.2$)')
plt.xlabel('t')
plt.ylabel('$\\upsilon$')
plt.plot(T, h)
plt.plot(T, ones,'r--', linewidth=0.8)
plt.plot(T, zeros, 'k--', linewidth = 0.8)
plt.savefig('nomedoarquivo.png')

#plot das animações

f = np.zeros(L, dtype = np.float32)
w = np.zeros(L, dtype = np.float32)
f[1:11] = 0.2
x = np.arange(L, dtype = np.float32)
t = 0
n = 1

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
frames = []
plt.title('Model FN, $\\ell$ = 10, I = 0.2')
plt.ylabel('$\\upsilon$')
plt.xlabel('x')
for m in range(1, 10000):

if (t == n):
curve = ax.plot(x, f, 'c')
frames.append(curve)
n = n + 1
aux = f[:]
f = Nagumo_CN(f, w, A)
w = FitzHung_FTCS(aux , w)

t = m*dt

animacao = ArtistAnimation(fig, frames, interval=50, blit=True)
animacao.save('nomedoarquivo.mp4')

</source>

===FTCS Modelo FitzHung-Nagumo 2D===
<source lang="python">

import numpy as np
import matplotlib
from matplotlib import pyplot as plt
import seaborn as sb

#Definindo Polinômio de Nagumo
def Pol(u):
a = 0.1
return u*(1 - u)*(u - a)

#Declarando Constantes
L = 100
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)

#Definindo função que calcula a equação de recuperação 2D
def FitzHung2D_FTCS(v, w):
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)
L = len(v)

aux = w[:]
for i in range(1, L-1):
for j in range(1, L-1):
aux[i,j] = w[i,j] + dt*0.005*(v[i,j] - 0.0001*w[i,j])

for i in range(1, L-1):
for j in range(1, L-1):
w[i,j] = aux[i,j]

return (w)

#Definindo função que calcula a equação de Nagumo em 2D como condições periódicas em y
def Nagumo2D_FTCS(v, w):
D = 0.1
dt = 1
dx = 1
r = (D*dt)/(dx**2)
L = len(v)

aux = v[:,:]
for j in range(1,L-1):
aux[0,j] = v[0,j] + r*(v[ L-1 ,j] + v[1,j] + v[0,j-1] + v[0,j+1] -4*v[0,j]) + dt*Pol(v[0,j]) - dt*w[0,j]

aux[0,L-1] = v[0,L-1] + r*(v[ L-1 ,L-1] + v[1,L-1] + v[0,L-2] + 0 -4*v[0,L-1]) + dt*Pol(v[0,L-1]) - dt*w[0,L-1]

for i in range(1, L-1):
for j in range(1, L-1):
aux[i,j] = v[i,j] + r*(v[i-1,j] + v[i+1,j] + v[i,j-1] + v[i,j+1] -4*v[i,j]) + dt*Pol(v[i,j]) - dt*w[i,j]

for j in range(1,L-1):
aux[L-1,j] = v[L-1,j] + r*(v[L-2,j] + v[0,j] + v[L-1,j-1] + v[L-1,j+1] -4*v[L-1,j]) + dt*Pol(v[L-1,j]) - dt*w[L-1,j]

aux[L-1,L-1] = v[L-1,L-1] + r*(v[L-2,L-1] + v[0,L-1] + v[L-1,L-2] + 0 -4*v[L-1,L-1]) + dt*Pol(v[L-1,L-1]) - dt*w[L-1,L-1]

for i in range(1, L):
for j in range(1, L-1):
v[i,j] = aux[i,j]

return (v)

#Função criada para evitar o uso de estruturas booleanas dentro do loop, que torna a simulação mais demorada
def Para_Evitar_Ifs(g):
intervalo = 20
T = g[0]
v = g[1]
w = g[2]
for t in range(T,T + intervalo):
aux = v[:,:]
v = Nagumo2D_FTCS(v, w)
w = FitzHung2D_FTCS(aux , w)

return ([T+intervalo,v,w])

#Criando os espaço e o grid de fotos
x = np.arange(L, dtype = np.float32)
y = np.arange(L, dtype = np.float32)
T = 0
dl = 4
p1 = (dl/2) - 1
p2 = (dl/2) + 1
tamanho = 2*L/dl

x,y = np.meshgrid(x,y)

f, (ax1,ax2,ax3,ax4) = plt.subplots(1,4,sharey=True, figsize=(24,4))
f.suptitle('Pulso Lateral de Tamanho '+str(int(tamanho))+'')

#Declarando condições iniciais e condições de contorno não periódicas
v = np.zeros((L,L))
v[int(p1*tamanho/2)-6:int(p2*tamanho/2),1:5] = 0.2

w = np.zeros((L,L))

#Plotando os 4 primeiros mapas de calor
g = [T,v,w]
graph1 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax1)
graph1.set_title('t='+str(g[0])+'')
graph1.invert_yaxis()

g = Para_Evitar_Ifs(g)
graph2 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax2)
graph2.set_title('t='+str(g[0])+'')
graph2.invert_yaxis()

g = Para_Evitar_Ifs(g)
graph3 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax3)
graph3.set_title('t='+str(g[0])+'')
graph3.invert_yaxis()

g = Para_Evitar_Ifs(g)
graph4 = sb.heatmap(g[1],cbar=True,vmin = -0.3, vmax = 1,ax=ax4)
graph4.set_title('t='+str(g[0])+'')
graph4.invert_yaxis()

plt.savefig('Lateral_1.png')
plt.show()

</source>

==Referências==
<references/>