izh_single_neuron_tonic_rk.m

It can be used in simulating izhikevich neurons

function x = izh_single_neuron_tonic_RK (I_ext)
a=0.02; b=0.2;  c=-65;  d=6;
V=-70;  u=b*V;
VV=[];  uu=[];
tau = 0.01; tau1=.01; tspan = 0:tau:30;
T1=tspan(end)/10;
for t=tspan
    if (t>T1) 
        I_ext=14;
    else
        I_ext=0;
    end;
    %V = V + tau*(0.04*V^2+5*V+140-u+I_ext);
     k1 = tau*f1(V, u, I_ext);
     k2 = tau*f1(k1/2+V, u, I_ext);
     k3 = tau*f1(k2/2+V, u, I_ext);
     k4 = tau*f1(k3+V, u, I_ext);
     V = V + (k1 + 2*k2 + 2*k3 + k4)/6;
    
    
     k1 = tau1*f2(V, u);
     k2 = tau1*f2(V, u+ (k1)/2);
     k3 = tau1*f2(V, u+ (k2)/2);
     k4 = tau1*f2(V, u + k3 );
     u = u + (k1 + 2*k2 + 2*k3 + k4)/6;
     
     %u = u + tau1*a*(b*V-u);
    if V > 30
        VV(end+1)=30;
        V = c;
        u = u + d;
    else
        VV(end+1)=V;
    end;
    uu(end+1)=u;
end;
plot(tspan,VV,[0 T1 T1 max(tspan)],-100+[0 0 10 10]);
axis([0 max(tspan) -100 30])
%axis off;
title(' Tonic spiking_Single_Neuron_Izhikevich_RK');

function dvbydt = f1(v,u,I)
dvbydt = (0.04*v.^2+5*v+140-u+I);

function du2bydt = f2(v,u)
a2=0.06; 
b2=0.22; 
du2bydt = a2*(b2*v-u);