Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Views: 88
%auto typeset_mode(True, display=False)

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:

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

Esta equação admite 3 equilíbrios:

var('x c') solve(c*(x-x^3/3),x)
(x,c)\left(x, c\right)
[x=3,x=3,x=0]\left[x = -\sqrt{3}, x = \sqrt{3}, x = 0\right]
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=±3x=\pm \sqrt{3}) e um instável (x=0x=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, yy, de "recuperação", que neutralize a excitação e leve o sistema de volta para o equilíbrio de "repouso".

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

dydt=1c(x+aby)\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 aa e bb são positivos e tomam valores preferencialmente nas seguintes faixas: 12b3<a<11-\frac{2b}{3} \lt a \lt 1 e 0<b<10 \lt b \lt 1. A forma como cc aparece em ambas as equações serve para ajustarmos a intensidade da excitabilidade em relação à recuperação. Aumentando cc, 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:

dydt=0y=1bx+ab\frac{dy}{dt}=0 \Rightarrow y = \frac{1}{b}x + \frac{a}{b}

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

 

Note que o parâmetro cc não afeta nenhuma das nullclines.

var('y') c=10 j=0 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(1/b*x+a/b,(-2,2),color='green',ymin=-1.5,ymax=1.5, legend_label="Nuliclina de $y$") show(vf+xnull+ynull)
yy
sols=solve([c*(x-1/3*x^3-y+j),1/c*(x+a-b*y)],[x,y]) for sol in sols: show(sol)
[x=(0.599704017622+1.35238113202i)\displaystyle x = \left(0.599704017622 + 1.35238113202i\right), y=(1.62463002203+1.69047641503i)\displaystyle y = \left(1.62463002203 + 1.69047641503i\right)]
[x=(0.5997040176221.35238113202i)\displaystyle x = \left(0.599704017622 - 1.35238113202i\right), y=(1.624630022031.69047641503i)\displaystyle y = \left(1.62463002203 - 1.69047641503i\right)]
[x=(1.19940806511)\displaystyle x = \left(-1.19940806511\right), y=(0.624260036595)\displaystyle y = \left(-0.624260036595\right)]
pe = point((-1.19940806511,-0.624260036595), color='red',pointsize=50) show(vf+xnull+ynull+pe)
def fun(t,Y, params): x, y = Y c,j,a,b=params #j=0.6*(heaviside(t-20)-heaviside(t-50)) return[ c*(x-x**3/3-y+j), (x+a-b*y)/c ] jfun = lambda t,st,params: [[10*st[0]^2,-10],[0.1,-0.08]] jfun (0,[-1,-.5],[.9,0,.7,.8])
[[10,10],[0.100000000000000,0.0800000000000000]]\left[\left[10, -10\right], \left[0.100000000000000, -0.0800000000000000\right]\right]
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.3,.7,.8],num_points=500)
#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
True\mathrm{True}
traj = list_plot([(p[1][0],p[1][1]) for p in T.solution],color="orange", plotjoined=True); show(vf+xnull+ynull+pe+traj)
%cython cimport sage.gsl.ode import sage.gsl.ode include 'gsl.pxi' cdef class nagumo_C(sage.gsl.ode.ode_system): cdef double j def __init__(self,j): self.j = j cdef int c_f(self,double t, double *z,double *dzdt): cdef: double x = z[0] double y = z[1] double c = 3 double j = self.j double a = 0.7 double b = 0.8 dzdt[0] = c*(x-x**3/3-y+j) dzdt[1] = (x+a-b*y)/c return GSL_SUCCESS
/projects/dc3261de-e369-4a08-bd7e-89894490f3e5/.sagemathcloud/sage_server.py:920: DeprecationWarning: Importing tmp_dir from here is deprecated. If you need to use it, please import it directly from sage.misc.temporary_file See http://trac.sagemath.org/17460 for details. code = code_decorator(code)
T = ode_solver() T.algorithm='rk8pd' cfun = nagumo_C(0.4) T.function= cfun t_range = [0,200] y0 = [-1.19940806511,-0.624260036595] T.ode_solve(t_range,y0,num_points=500)
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);
import numpy as np @interact def bifurcation(j=slider(0,1,0.01,0)): T = ode_solver() T.algorithm='rk8pd' pts = [] cfun = nagumo_C(j) T.function= cfun t_range = [0,100] y0 = [-1.19940806511,-0.624260036595] T.ode_solve(t_range,y0,num_points=1000) sol = T.solution 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, gridlines=True);
Interact: please open in CoCalc
import numpy as np @interact def bifurcation(j_min=slider(0,1,0.01,0.3),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,500): cfun = nagumo_C(j) T.function= cfun t_range = [0,50] y0 = [-1.19940806511,-0.624260036595] T.ode_solve(t_range,y0,num_points=500) sol = T.solution pts += [(j,p[1][0]) for n,p in enumerate(sol[0:]) if p[1][1]-0.5 < 0.01] show(points(pts,pointsize=2), gridlines=True)
Interact: please open in CoCalc