% fitzhugh_nagumo.m

%  NEUR 1680
%  Last updated March 7, 2012

% V' = V - V^3/3 - W + I
% W' = 0.08*(V + 0.7 - 0.8*W)

% This script allows you to do a (T,D) analysis:
%   after you complete the trajectory-creation step,
%   click at the equilibrium point, to let the code calculate
%   and display the values of T and D at that point.


clear

min_V = -2;
max_V = 2;
V_step = .0001;
min_W = -1.5;
max_W = 2.5;
dt = .001;      % time step
t_max = 100;     % length of simulation
line_width = 2;

I = .36;    % Start at I = -0.3 and increase to 1.5
            % Examine I = 0
            % Note Hopf bifurcations around 0.34 and around 1.42

            
% Use the following parameters to perform a detailed analysis
%  near the bifurcation (e.g. I = 0.33 and I = 0.35)
% min_V = -1.5;
% max_V = -0.5;
% V_step = .0001;
% min_W = -0.5;
% max_W = 0;
% dt = .001;      % time step
% t_max = 1000;     % length of simulation
% line_width = 1;


VV = min_V:V_step:max_V;    % V axis

figure(1)

% NULLCLINES
Vertical_Nullcline = VV - VV.^3/3 + I;  
Horizontal_Nullcline = (VV + 0.7)/0.8; 
plot(VV,Vertical_Nullcline,'g',VV,Horizontal_Nullcline,'m')      
axis([min_V max_V min_W max_W])
grid on
hold on

% DIRECTION FIELD 
% Define the grid points:
[x,y] = meshgrid(linspace(min_V,max_V,15),linspace(min_W,max_W,15));
F = @(x,y) x - x.^3/3 - y + I ;
G = @(x,y) 0.08.*(x + 0.7 - 0.8.*y) ;

% The derivatives at all grid points:
dx = F(x,y);
dy = G(x,y);

% Put the (dx,dy) arrow at each grid point (x,y):
quiver(x,y,dx,dy)



% To create a new trajectory, click mouse on desired initial point in (x,y) plane.

% To end simulation (e.g. to change parameter values),
%   click mouse in grey margin, i.e., outside of permissible area.

time_axis = 0:dt:t_max;
n_time_steps = length(time_axis);

V = .5*(min_V + max_V)*ones(1,n_time_steps);
W = .5*(min_W + max_W)*ones(1,n_time_steps);
    % will contain an entire trajectory

while V(1)>min_V && V(1)<max_V && W(1)>min_W && W(1)<max_W  
    % will create a new trajectory until mouse is clicked in margin

    figure(1)
    [V(1),W(1)] = ginput(1);    % click mouse on desired initial value

    for i = 2:n_time_steps
        V(i) = V(i-1) + (V(i-1) - V(i-1)^3/3 - W(i-1) + I)*dt;
        W(i) = W(i-1) + (0.08*(V(i-1) + 0.7 - 0.8*W(i-1)))*dt;
    end
    plot(V,W,'r', 'LineWidth',line_width)
    hold on
    drawnow

    figure(2)
    plot(time_axis,V,time_axis,W)
    grid on
    xlabel('time')
    legend('membrane potential', 'recovery variable')
    title(['INPUT CURRENT: ',num2str(I)])

end

hold off


% T-D ANALYSIS OF LINEARIZED SYSTEM:

figure(1)
title(['INPUT CURRENT: ',num2str(I)])
xlabel('membrane potential V')
ylabel('recovery variable W')
[V0,W0] = ginput(1);    % click mouse at equilibrium point
hold off

T = .936 - V0^2;
D = .016 + .064*V0^2;

% Visualize the parameter set in the T-D plane:

figure(3)
plot(T,D,'*r','MarkerSize',20)

% We also plot the two axes, and the parabola of equation D = .25*T^2:

hold on
min_T = -1;
max_T =  1;
min_D = -.2;
max_D =  .2;
plot([min_T max_T],[0 0])
plot([0 0],[min_D max_D])
T_axis = min_T:.01:max_T;
plot(T_axis,.25*T_axis.^2)
axis([min_T max_T min_D max_D])
grid on
hold off

xlabel('TRACE')
ylabel('DETERMINANT')
title(['INPUT CURRENT    I = ',num2str(I)])
