%% HW4 - The Hodgkin and Huxley model - FRTF01
% Department of Automatic Control, Lund University

%% Section 1 - Constants
% Na - sodium, K - potassium, L - Leakage

% Membrane capacitance [microF/cm^2]
C_m  = % fill in
% Maximum conductances [mS/cm^2]
g_Na =  % fill in
g_K =   % fill in
g_L  = % fill in
% Nernst reversal potentials [mV]
E_Na = % fill in
E_K  = % fill in
E_L  = % fill in

%% Section 2 - Channel gating kinetics, rate functions
% Functions of membrane voltage
alpha_m = @(V) 0.1*(V+45)./(1-exp(-(V+45)./10));
beta_m = @(V) 4*exp(-(V+70)./18);
alpha_h = @(V) 0.07*exp(-(V+70)./20);
beta_h = @(V) 1./(1+exp(-(V+40)./10));
alpha_n = @(V) 0.01*(V+60)./(1-exp(-(V+60)./10));
beta_n = @(V) 0.125*exp(-(V+70)./80);

%% Section 3 - Membrane currents (in microA/cm^2)
I_L = @(V) % fill in
I_Na = @(V,m,h) % fill in
I_K = @(V,n) % fill in

%% Section 4 - External current
I_ext = @(t) % fill in


%% Section 5 - Defining the differential equation
% Define X = [V,m,n,h];

dAlldt = @(t,X) ...
    [(I_ext(t)-I_Na(X(1),X(2),X(4))-I_K(X(1),X(3))-I_L(X(1)))./C_m;...
    alpha_m(X(1))*(1-X(2)) - beta_m(X(1))*X(2);...
    alpha_n(X(1))*(1-X(3)) - beta_n(X(1))*X(3);...
    alpha_h(X(1))*(1-X(4)) - beta_h(X(1))*X(4)];

%% Section 6 - The steady state

% Initial values = [V(0), m(0), n(0), h(0)]
initial_values = [0 0 0 0]; % Just starts somewhere to find out the steady state.
tspan = [0 500]; % Calculates the solution during 500 [ms]

[t, X] = ode45(dAlldt, tspan, initial_values);

V = X(:,1); % 1st column is the V values
m = X(:,2); % 2nd column is the m values
n = X(:,3); % 3rd column is the n values
h = X(:,4); % 4th column is the h value

figure(1)
subplot(3,1,1);
plot(t, V, 'k');
title('1A: Find steady state values')
ylabel('V (mV)');
xlabel('t (ms)')

subplot(3,1,2);
plot(t, m, 'r');
hold on
plot(t, n, 'g');
hold on
plot(t, h, 'b');
ylabel('Gating Value')
xlabel('t (ms)')
legend('m', 'n', 'h', 'Location','SouthEast')

subplot(3,1,3);
plot(t, I_ext(t), 'k');
xlabel('t (ms)');
ylabel('I_{ext} (\mu{A}/cm^2)');
ylim([-1 max(I_ext(t))+1]);


%% Section 7 -  Find the threshold potential and specify the refractory period

initial_values = []; % Fill in the ones you determined in the section 6
k = % use different k-values to determine the threshold potential
I_ext = @(t) k.*floor(t./100); % use help floor to understand what it does

% Run the section to create the figure with the new I_ext

dAlldt = @(t,X) ...
    [(I_ext(t)-I_Na(X(1),X(2),X(4))-I_K(X(1),X(3))-I_L(X(1)))./C_m;...
    alpha_m(X(1))*(1-X(2)) - beta_m(X(1))*X(2);...
    alpha_n(X(1))*(1-X(3)) - beta_n(X(1))*X(3);...
    alpha_h(X(1))*(1-X(4)) - beta_h(X(1))*X(4)];

[t, X] = ode45(dAlldt, tspan, initial_values);

V = X(:,1); % the first column is the V values
m = X(:,2); % the second column is the m values
n = X(:,3); % the third column is the n values
h = X(:,4); % the fourth column is the h value

figure(2)
subplot(3,1,1);
plot(t, V, 'k');
title('Find the threshold potential')
ylabel('V (mV)');
xlabel('t (ms)')

subplot(3,1,2);
plot(t, m, 'r');
hold on
plot(t, n, 'g');
hold on
plot(t, h, 'b');
ylabel('Gating Value')
xlabel('t (ms)')
legend('m', 'n', 'h', 'Location','SouthEast')

subplot(3,1,3);
plot(t, I_ext(t), 'k');
xlabel('t (ms)');
ylabel('I_{ext} (\mu{A}/cm^2)');
ylim([-1 max(I_ext(t))+1]);
