%%  Part 1:  Download and examine the HH model.
%   To start, download the function:
%
%   HH.m
%
%   from http://math.bu.edu/people/mak/samsi/HH.m
%
%   Examine this function ... is it consistent with the model we proposed? 

%%  Part 2:  At low I, examine the HH dynamics.  %%%%%%%%
%   To understand how the HH model works, we'll start by focusing on the
%   case when I is small. Let's fix the input current to zero,

I = 0;

%   and simulate the model for 100 ms,

T = 100;

%   We've now defined both inputs to the "HH" function, and can execute it,
%   as follows,

[V,m,h,n,t] = HH(I,T);

%   Notice that the function returns five outputs, which we assign to the
%   variables 'V', 'm', 'h', 'n', and 't'.

%Q:  What are the dynamics of the model variables resulting from this
%simulation?

%A:  Let's make some plots.

subplot(2,1,1)                  %Make a subplot,
plot(t,V,'k')                   %... and plot the voltage,
ylabel('V [mV]')                %... with y-axis labeled.

%  Then, to plot the gating variables, 

subplot(2,1,2)                  %Make a second subplot,
plot(t,m,'r')                   %... and plot the sodium activation variable,
hold on                         %... hold the graphics window,
plot(t,h, 'b')                  %... and plot the sodium inactivation variable,
plot(t,n, 'g')                  %... and plot the potassium activation variable .
hold off                        %... release the graphics window,
xlabel('Time [ms]')             %... label the x-axis.
ylabel('mS/cm^2')               %... and label the y-axis.
legend('m', 'h', 'n')

%ASIDE:  Here's a nice plotting trick, to link the x-axes of our two
%subfigures.  Linking the axes is useful so that, when we zoom or move one
%subfigure, the other subfigure will match the x-axis.

ax(1) = subplot(2,1,1);
ax(2) = subplot(2,1,2);
linkaxes([ax(2) ax(1)], 'x')

%CONCLUSION:  At this value of input current (I=0), the model dynamics
%approach a "fixed point", whose location we can identify as a point in 
%four dimensional space.

%%  Part 3:  At high I, examine the HH dynamics of a spike.  %%%%%%%%%%
%   Let's now increase the input current to the HH model and get this model
%   to generate repeated spiking activity.  To do so, let's set,

I = 10;

%   We can now simulate this model,

[V,m,h,n,t] = HH(I,T);

%IN LAB Q:  What happens to the dynamics? 

%IN LAB A:  Let's make some plots.

subplot(2,1,1)                  %Make a subplot,
plot(t,V,'k')                   %... and plot the voltage,
xlim([42 56])                   %... focus on a single spike,
ylabel('V [mV]')                %... with y-axis labeled.

%  Then, to plot the gating variables, 

subplot(2,1,2)                  %Make a second subplot,
plot(t,m,'r')                   %... and plot the sodium activation variable,
hold on                         %... hold the graphics window,
plot(t,h, 'b')                  %... and plot the sodium inactivation variable,
plot(t,n, 'g')                  %... and plot the potassium activation variable .
hold off                        %... release the graphics window,
xlim([42 56])                   %... focus on a single spike,
xlabel('Time [ms]')             %... label the x-axis.
ylabel('mS/cm^2')               %... and label the y-axis.
legend('m', 'h', 'n')

ax(1) = subplot(2,1,1);
ax(2) = subplot(2,1,2);
linkaxes([ax(2) ax(1)], 'x')

%IN LAB Q:  How do the conductances evolve during a spike?

%% CHALLENGES:
%
%  1. Compute and display the f-I curve for the Hodgkin-Huxley model.
%  2. How do the model dynamics change when the maximal sodium conductance is
%  reduced, or enhanced?
%  3. How do the model dynamics change when the maximal potassium conductance is
%  reduced, or enhanced?