EDO - NEURO - SEANCE 2 - Modeles Integre-et-Tire et Modeles de Phase

Contents

Modele d'Izhikevich, methode d'Euler

Les equations:

C*dv/dt = k*( v - vr )*( v - vt ) - u + I
du/dt = a*( b*( v - vr ) - u )

Si v>=vpeak alors v=c, u=u+d

% parametres du modele
C=100; vr=-60; vt=-40; k=0.7;
a=0.03; b=-2; c=-50; d=100;
vpeak=35;

% parametres d'integration
T=1000; dt=1;
n=round(T/dt);
v=vr*ones(1,n); u=0*v;
I=[zeros(1,0.1*n),70*ones(1,0.9*n)];

% boucle principale: methode d'Euler
for i=1:n-1
    % Integre (etape Euler)
    v(i+1) = v(i) + dt * (k*(v(i)-vr)*(v(i)-vt)-u(i)+I(i))/C;
    u(i+1) = u(i) + dt * a*(b*(v(i)-vr)-u(i));

    % Tire (reinitialisation)
    if v(i+1) >= vpeak
        v(i)   = vpeak;
        v(i+1) = c;
        u(i+1) = u(i+1)+d;
    end
end

plot(dt*(1:n),v,'g');

Modele d'Izhikevich, methode 'Event Locator'

La fonction izhikevich implemente le modele integre-et-tire d'Izhikevich.

sol = izhikevich;

Le code de la fonction izhikevich est

function sol = izhikevich
% IZHIKEVICH modele neurone integre-et-tire

% parametres des equations
C=100;
vr=-60;
vt=-40;
k=0.7;
a=0.03;
b=-2;

intensite = 70;

vpeak=35;
c=-50;
d=100;

% options de ode
% l'option 'Events' permet de gerer des evenements
% les ... permettent de continuer l'expression sur la ligne suivante
options = odeset('Events',@it_event, ...
               'AbsTol',1e-6, ...
               'RelTol',1e-3);

% conditions initiales
CI = [vr; 0];

% parametres de simulation
tfinal = 1000;
t0 = 0;
tout = t0;
xout = CI';
teout = [];
xeout = [];

while t0<tfinal
    [t,x,te,xe] = ode45(@it,[t0 tfinal],CI,options);
    nt = length(t);

     % on accumule les solutions
     tout = [tout; t(2:nt)];
     xout = [xout; x(2:nt,:)];
     teout = [teout; te];
     xeout = [xeout; xe];
     % on met a jour t0
     t0 = t(nt);
     % on met a jour les CI
     CI = [c; x(nt,2)+d];
 end
% on affiche la solution pour v [x(1)]
plot(tout,xout(:,1),'r')

% on met la solution dans une structure de sortie
sol.t = tout;
sol.x = xout;
sol.te = teout;
sol.xe = xeout;

%--------------------------------------------------------
% Fonctions imbriquees
%--------------------------------------------------------

      function dxdt = it(t,x)
      % IZ equations EDO
         v=x(1);
         u=x(2);
         dxdt = [(k*(v-vr).*(v-vt)-u+input_current(t))/C;
               a*(b*(v-vr)-u)];
     end
     function I = input_current(t)
     % INPUT_CURRENT courant applique a la membrane
         I = intensite * ( t>=(0.1*tfinal) );
     end
     function [value,isterminal,direction] = it_event(~,x)
     % IZ_EVENT detecte les evenements
         value = x(1)-vpeak;
         isterminal = 1;
         direction = 1;
     end
end

Modele de phase

On peut identifier un point sur un cycle limite avec une phase, entre 0 et $2 \pi$.

Le modele de phase le plus simple, represente en coordonnees polaires est

dphi/dt = f(phi);

si phi>=2*pi, alors phi = 0

avec f(phi) (strictement) positive.

La phase est discontinue, mais, representee en coordonnees cartesiennes, on obtient une hauteur y = sin(phi) continue. Pour completer, on peut tracer le portrait de phase avec x = cos(phi) et y = sin(phi).

Modele de phase avec relaxation

En plus de la phase, on peut inclure l'amplitude comme variable dynamique. Cette deuxieme variable, r, est positive et est regie par l'equation

dr/dt = g(r,p)

L'equation de phase devient dependante de l'amplitude

dp/dt = f(r,p);

si p>=2*pi, alors p = 0

On obtient alors, en coordonnees cartesiennes, x = r*cos(p) et y = r*sin(p).

Exercice 1

Comparez les solutions du modele d'Izhikevich obtenues avec la methode d'Euler et la methode 'Event Locator'. Quelle methode est la plus precise? Comment s'assurer que la solution obtenue est precise? Comparez les temps d'execution de chaque methode. Variez le pas de temps dt et les les parametres AbsTol et RelTol.

Exercice 2

En utilisant la methode de resolution de votre choix, tracez le portrait de phase du modele d'Izhikevich (tracez u contre v). Localisez les points de discontinuite.

Exercice 3

Codez une methode pour resoudre l'equation de phase suivante:

dp/dt = w*(1.01+sin(phi))

Utilisez le parametre w=1.0, la condition initiale CI = 0.0, pour des temps entre 0 et 200. Pour la methode numerique, inspirez-vous, c.-a-d. copiez et modifiez, une des deux methode pour le modele d'Izhikevich. Tracer la solution en coordonnees cartesiennes. Utilisez la fonction comet pour tracer la trajectoire dans le plan de phase (x,y).

Exercice 4

a) Codez une methode pour resoudre l'equation de phase suivante:

dr/dt = - lambda*(r-A)

dphi/dt = omega

avec les parametres lambda = 0.001, 0.01, 0.1, A = 1, et omega = 2*pi/100. entre t = 0 et t = 1000. Resolvez pour des conditions initiales phi0 = 0, r0 = 0.1, 1.0, 2.0. Tracer les solutions en coordonnees cartesiennes.

b) Modifiez les equations pour inclure un terme de forcage dans la direction des x en coordonnees cartesiennes:

dr/dt = - lambda*(r-A) + F*cos(phi)*cos(2*pi/T*t)

dphi/dt = omega - 1/r*F*sin(phi)*cos(2*pi/T*t)

avec T = 50, 100, 200, et F assez petit.


Published with MATLAB® R2012b