%FHN_3D_FIGURE.M: Evolution of Fitzhugh-Nagumo equations plotted in 3D Space
%Sitabhra Sinha

clear all
%PARAMETERS 
%[from Glass & Josephson PRL 75 (1995) 2059]
%alpha=0.139;gamma=0.394;silonn=0.005; %0.02;
%from Singh et al, PRL (2012)
A=5.0;alpha=0.2;silonn=0.05;gamma=1;
%External current
%Iext=0.1; %+ve: Depolarizing current
Iext=-0.5; %-ve: Hyperpolarizing current
%INITIAL CONDITIONS
u=0.0;w=0.0;
%INTEGRATION PARAMETER
dt=0.01;

%Stimulation parameters
t_on=4000;t_off=5500; %4500;

%TIME EVOLUTION
for i=1:20000,
%du=u*(1-u)*(u-alpha)-w+0.1;
du=A*u*(1-u)*(u-alpha)-w+Iext*(0.5+0.5*sign(i-t_on))*(0.5+0.5*sign(t_off-i));
dw=silonn*(gamma*u-w); %-0.075); %-0.075);
urec(i)=u+dt*du;wrec(i)=w+dt*dw;
u=urec(i);w=wrec(i);
end;

figure, plot3([1:length(urec)]*dt,urec,wrec,'k'),
xx=[-0.4:0.01:1];hold on, plot3(zeros(size(xx)),xx,A*xx.*(1-xx).*(xx-alpha),'r :')
xx=[-0.2:0.01:0.8];plot3(zeros(size(xx)),xx,gamma*xx,'b :')
plot3(zeros(size(urec)),urec,wrec,'k')
plot3([1:length(urec)]*dt,urec,0.8*ones(size(wrec)),'r')
plot3([1:length(urec)]*dt,ones(size(urec)),wrec,'b')
%stimulation
plot3([1:length(urec)]*dt,Iext*(0.5+0.5*sign([1:length(urec)]-t_on)).*(0.5+0.5*sign(t_off-[1:length(urec)])),-0.05*ones(size(wrec)),'r')

axis([0 i*dt -0.5 1.2 -0.2 0.8])
view([20 22])
set(gca,'FontSize',14)
xlabel('Time (a.u.)')
ylabel('Activation u')
zlabel('Inactivation v')
set(get(gca,'xlabel'),'rotation',-3);
set(get(gca,'ylabel'),'rotation',40);
box on