Aulas do Curso de Modelagem matemática IV da FGV-EMAp

Path: Modelagem-Matematica-IV/Planilhas Sage / Aula 8 - O Modelo de Fitzhugh-Nagumo.ipynb

Views: ²⁴¹²
License: GPL3
Image: default

Kernel: SageMath 9.3

Modelo de Fitzhugh-Nagumo

In [16]:

%display typeset

O modelo de Fitzhugh Nagumo é um elegante modelo de um sistema excitável. Históricamente foi desenvolvido de um esforço por simplificar o modelo original de Hodgkin-Huxley mantendo suas propriedades dinâmicas intactas. Mas este modelo é um ponto de partida para o estudo do fenômeno da excitabilidade que, sendo bastante comum em biologia, não se restringe à eletrofisiologia dos neurônios.

Computador analógico usado por Richard Fitzhugh.

Richard Fitzhugh com o computador analógico usado para calcular o modelo.

O que é excitabilidade?

Um sistema excitável pode ser descrito como um sistema dinâmico que possui um estado de repouso para o qual sempre retorna após sofrer pequenas pertubações. Entretanto se a perturbação ou estímulo ultrapassar seu limiar de excitabilidade, o sistma fará uma excursão mais longa pelo espaço de estados antes de retornar ao seu estado de repouso. Durante esta escursão será refratário a novos estímulos, ou seja, novos estímulos não afetarão sua trajetória. Mas após retornar ao repouso, estará sujeito a uma nova excitação.

O modelo de Fitzhugh-Nagumo consiste em duas equações apenas. A primeira busca representar a excitação do sistema:

$\frac{dx}{dt} = c \left(x-\frac{1}{3}x^3 \right)$

Esta equação admite 3 equilíbrios:

In [1]:

var('c')
solve(c*(x-x^3/3),x)

[x == -sqrt(3), x == sqrt(3), x == 0]

In [2]:

p=plot(1*(x-x^3/3),(-2,2))
po = point([(-sqrt(3),0),(0,0),(sqrt(3),0)],color='red',pointsize=40)
show(p+po)

Pelo gráfico, podemos ver que temos dois equilíbrios estáveis ( $x=\pm \sqrt{3}$ ) e um instável ( $x=0$ ). Mas esta equação sozinha, nos dá uma sistema bistável, e não um excitável. Para isso temos que acrescentar uma variável, $y$ , de "recuperação", que neutralize a excitação e leve o sistema de volta para o equilíbrio de "repouso".

$\frac{dx}{dt} = c \left(x-\frac{1}{3}x^3 -y + j \right)$

$\frac{dy}{dt}=\frac{1}{c}(x+a-by)$

o parâmetro j representa o estímulo externo, no caso do neurônio a correte de despolarização injetada na célula. Os parâmetros $a$ e $b$ são positivos e tomam valores preferencialmente nas seguintes faixas: $1-\frac{2b}{3} \lt a \lt 1$ e $0 \lt b \lt 1$ . A forma como $c$ aparece em ambas as equações serve para ajustarmos a intensidade da excitabilidade em relação à recuperação. Aumentando $c$ , aumentamos a excitabilidade e diminuimos a recuperação.

Análise no Plano de Fase

Para nos ajudar nesta análise, vamos definir as nuliclinas de ambas as variáveis:

$\frac{dy}{dt}=0 \Rightarrow y = \frac{1}{b}x + \frac{a}{b}$

$\frac{dx}{dt}=0 \Rightarrow y= x- \frac{1}{3}x^3 +j$

Note que o parâmetro $c$ não afeta nenhuma das nuliclinas.

In [10]:

var('x y a b c j')
dxdt = c*((x-(x**3)/3)-y+j)
solve(dxdt,y)

[y == -1/3*x^3 + j + x]

In [11]:

dydt = (x+a-b*y)/c
solve(dydt,y)

[y == (a + x)/b]

In [12]:

var('y')
c=10
j=0.3
a=.7
b=.8
vf = plot_vector_field([c*((x-(x**3)/3)-y+j),(x+a-b*y)/c],(x,-2,2),(y,-1.5,1.5),axes_labels=[r'$x$',r'$y$'])
xnull = plot(x-x^3/3+j,(-2,2),color='blue',ymin=-1.5,ymax=1.5, legend_label="Nuliclina de $x$")
ynull = plot((a+x)/b,(-2,2),color='green',ymin=-1.5,ymax=1.5, legend_label="Nuliclina de $y$")
show(vf+xnull+ynull)

A interseção entre as nuliclinas é um equilíbrio do systema. Podemos calcular o seu valor.

In [13]:

sols=solve([dxdt(a=a,b=b,c=c,j=j),dydt(a=a,b=b,c=c)],[x,y])
for sol in sols:
    show(sol)

\renewcommand{\Bold}[1]{\mathbf{#1}}\left[x = \left(0.4966487372747987 + 1.220647330194089i\right), y = \left(1.495810921593498 + 1.525809162742611i\right)\right]

\renewcommand{\Bold}[1]{\mathbf{#1}}\left[x = \left(0.4966487372747987 - 1.220647330194089i\right), y = \left(1.495810921593498 - 1.525809162742611i\right)\right]

\renewcommand{\Bold}[1]{\mathbf{#1}}\left[x = \left(-0.9932974689126025\right), y = \left(-0.3666218551949299\right)\right]

In [14]:

pe = point((sol[0].rhs(),sol[1].rhs()), color='red',pointsize=50)
show(vf+xnull+ynull+pe)

Vamos calcular a Jacobiana do sistema:

In [17]:

Jac = jacobian([dxdt,dydt], [x,y])
Jac

\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} -{\left(x^{2} - 1\right)} c & -c \\ \frac{1}{c} & -\frac{b}{c} \end{array}\right)

In [18]:

eqx = -0.9932974689126025
Jac(x=eqx,b=b,c=c).eigenvalues()

\renewcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{1470389350} i \, \sqrt{2137383719017539699} + \frac{39407451}{1470389350}, \frac{1}{1470389350} i \, \sqrt{2137383719017539699} + \frac{39407451}{1470389350}\right]

In [19]:

def fun(t,Y, params):
    x, y = Y
    c,j,a,b=params
    #j=0.9*(heaviside(t-20)-heaviside(t-50))
    return[
    c*(x-x**3/3-y+j),
    (x+a-b*y)/c
    ]
jfun = lambda t,st,p: jacobian([c*((x-(x**3)/3)-y+j),(x+a-b*y)/c], [x,y]).subs(x=st[0], c=p[0],j=p[1],a=p[2],b=p[3])
jfun (0,[sol[0].rhs(),sol[1].rhs()],[c,j,a,b])

\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 0.1336013825181741 & -10 \\ \frac{1}{10} & -0.0800000000000000 \end{array}\right)

In [30]:

T = ode_solver()
T.algorithm='rk8pd'
T.function= fun
T.jacobian = jfun
t_range = [0,80]
y0 = [-1.19940806511,-0.624260036595]#[.5,.5]
T.ode_solve(t_range,y0,params=[3,0.5,.7,.8],num_points=500)

In [31]:

#fun(0,[1.49900000000000,0.0833333333333333],[0.08,0,.7,.8])
px = list_plot([(p[0],p[1][0]) for p in T.solution], plotjoined=True, legend_label="x");
py = list_plot([(p[0],p[1][1]) for p in T.solution], color='green', legend_label="y");
px.legend()
show(px+py);
#print T.solution

In [32]:

@interact
def traj_plot(x0=slider(-2,2,.1, -1.19940806511),y0=slider(-1.5,1.5,.1,-0.624260036595)):
    inits = [x0,y0]#[.5,.5]
    T.ode_solve(t_range,inits,params=[c,j,a,b],num_points=500)
    traj = list_plot([(p[1][0],p[1][1]) for p in T.solution],color="orange", plotjoined=True);
    pi = point([inits[0],inits[1]],color='green', pointsize=50)
    show(vf+xnull+ynull+traj+pi)

In [33]:

T = ode_solver()
T.algorithm='rk8pd'
#cfun = nagumo_C(0.4)
T.function= fun
t_range = [0,80]
y0 = [-1.19940806511,-0.624260036595]
T.ode_solve(t_range,y0,params=[3,0.4,.7,.8],num_points=500)

In [34]:

px = list_plot([(p[0],p[1][0]) for p in T.solution], plotjoined=True, legend_label="x");
py = list_plot([(p[0],p[1][1]) for p in T.solution], color='green', legend_label="y");
show(px+py);

In [35]:

import numpy as np
@interact
def bifurcation(j_min=slider(0,1,0.01,0),j_max=slider(0,2,0.01,1.4)):
    T = ode_solver()
    T.algorithm='rk8pd'
    pts = []
    for j in np.linspace(j_min,j_max,200):
        #cfun = nagumo_C(j)
        T.function= fun
        t_range = [0,50]
        y0 = [-1.19940806511,-0.624260036595]
        T.ode_solve(t_range,y0,params=[3,j,.7,.8], num_points=200)
        sol = T.solution
        pts += [(j,p[1][0]) for n,p in enumerate(sol[0:]) if abs(p[1][1]-0.5) < 0.01]
    show(points(pts,pointsize=2), gridlines=True)

In [ ]:

In [ ]: