% This is an implementation of the basic HH equations. It allows the user
% to provide a single stimulus at the time, strength, and duration of the
% user's choosing.

% Equations from Hodgkin, Huxley, J. Physiol (1952) 117, 500-544

% Iyad Obeid
% 25 November 2015
% Version 0.0
clc;
clear;
clf;

%%%%%%%%%%%%%%%%%%
% DEFINE CONSTANTS
%%%%%%%%%%%%%%%%%%




% Nernst Potentials
Ena = 52.4; Ek = -72.1; El = -49.2;

% Maximum Conductances
gna = 120; gk = 36; gl = 0.3;

% Membrane Capacitance
C = 1;

Re = 25;
Ri = 35;



a=0.01;
b=20;
r=@(z)a*z+b
drdz = a;
dsdz=1/sqrt(1+drdz^2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINE VOLTAGE-DEPENDENT GATE ACTIVATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% % Likelihoods of gates opening
% an = @(u) (0.1-0.01*u)/(exp(1-0.1*u)-1);
% am = @(u) (2.5-0.1*u)/(exp(2.5-0.1*u)-1);
% ah = @(u) 0.07*exp(-u/20);
% 
% % Likelihoods of gates closing
% bn = @(u) 0.125*exp(-u/80);
% bm = @(u) 4*exp(-u/18);
% bh = @(u) 1/(exp(3-0.1*u)+1);

% Likelihoods of gates opening
am = @(u) 0.1*(u+40)/(1-exp(-0.1*(u+40)));
an = @(u) 0.01*(u+55)/(1-exp(-0.1*(u+55)));
ah = @(u) 0.07*exp(-0.05*(u+65));

% Likelihoods of gates closing
bm = @(u) 4*exp(-0.0556*(u+65));
bn = @(u) 0.125*exp(-0.0125*(u+65));
bh = @(u) 1/(1+exp(-0.1*(u+35)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINE FORMULAE FOR STEADY STATE GATE ACTIVATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m_inf = @(u) am(u) / ( am(u) + bm(u) );
n_inf = @(u) an(u) / ( an(u) + bn(u) );
h_inf = @(u) ah(u) / ( ah(u) + bh(u) );

%%%%%%%%%%%%%
% DEFINE TIME
%%%%%%%%%%%%%

dt = 0.01;
T=10;
t = 0:dt:T;
nt=length(t);

%%%%%%%%%%%%%
% DEFINE length
%%%%%%%%%%%%%

L=3;
N=11;
dx=L/N;
l=0:dx:L;
nl=length(l);

dt/dx^2

%%%%%%%%%%%%%
% DEFINE parametrs
%%%%%%%%%%%%%

alpha = 1/(C*dx^2);
betta = 1/(2*C*dx);
R=1/(2*(Re+Ri));

%%%%%%%%%%%%%
% DEFINE initial condition
%%%%%%%%%%%%%
vse(1:nt, 1:nl) = -60*ones;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINE STIMULUS STRENGTH, DURATION, & DELAY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tSTIM_START = 2;
tSTIM_DUR = 1;
STIM_STRENGTH = 1000;

tSTIM_START1 = 30;
tSTIM_DUR1 = 1;
STIM_STRENGTH1 = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INITIALIZE STATE VARIABLES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m = 0*t + m_inf(0);
n = 0*t + n_inf(0);
h = 0*t + h_inf(0);
v = 0*t;
v(1)=-56;



%%%%%%%%%%%%%
% DEFINE MISC
%%%%%%%%%%%%%

inRange = @(x,a,b) (x>=a) & (x<b);

%%%%%%%%%%%
% MAIN LOOP
%%%%%%%%%%%



for i = 1:length(t)-1
%     for j = 2:nl-1
    j=2;
    % extract membrane voltage
    u = v(i);

    % solve for membrane currents
    Ik  = gk  * n(i)^4      * (u - Ek);
    Ina = gna * m(i)^3*h(i) * (u - Ena);
    Il  = gl  *               (u - El);   
    Imem = Ik + Ina + Il;
    
    % determine stimulus current, if any
    Istim = 0;
    if inRange( t(i) , tSTIM_START , tSTIM_START+tSTIM_DUR)
        Istim = STIM_STRENGTH;
    end
    
    if inRange( t(i) , tSTIM_START1 , tSTIM_START1+tSTIM_DUR1)
        Istim = STIM_STRENGTH1;
    end

    % define the state variable derivatives
    dv = (Istim - Imem)/C;
    dm = am(u) * (1-m(i)) - bm(u) * m(i);
    dh = ah(u) * (1-h(i)) - bh(u) * h(i);
    dn = an(u) * (1-n(i)) - bn(u) * n(i);
    
    % use forward euler to increment the state variables
    v(i+1) = v(i) + dv*dt;
    m(i+1) = m(i) + dm*dt;
    h(i+1) = h(i) + dh*dt;
    n(i+1) = n(i) + dn*dt;
 
    
    vse(i,j-1)=v(i+1);
%     vse(i,nl)=v(i+1);
    

%     vse(i,1) = v(i);
%     vse(1,:) = v(i);
%     vse(i,nl) = v(i);
    

     for j = 2:nl-1
%         vse(i+1,j) = vse(i,j)+alpha*R*(vse(i,j-1)-2*vse(i,j)+vse(i,j+1));
%         vse(i+1,j) = vse(i,j)+dt*(alpha*R*(vse(i,j-1)-2*vse(i,j)+vse(i,j+1))+dv);


%         vse(i+1,j) = vse(i,j)+dt*r(j)*(alpha*R*(vse(i,j-1)-2*vse(i,j)+vse(i,j+1))); %Рабочее
 
%           vse(i+1,j) = vse(i,j)+dt*(r(j)*R*alpha*(vse(i,j-1)-2*vse(i,j)+vse(i,j+1))+drdz*2*R*betta*(vse(i,j-1)+vse(i,j+1))); %Рабочее
          
           vse(i+1,j) = vse(i,j)+dt*dsdz*(r(j)*R*alpha*(vse(i,j-1)-2*vse(i,j)+vse(i,j+1))+drdz*2*R*betta*(vse(i,j-1)+vse(i,j+1)));
        
    end
    
end


% vse(length(t),:) = [];
% t(length(t))=[];
% v(length(t))=[];

% plot the results
figure(1)
plot(t,v);
xlabel('time (ms)');
ylabel('membrane voltage (mV)');
grid on

figure(2);
s=surf(vse);
xlabel('Length');
ylabel('Time');
zlabel('Voltage');
colorbar
s.EdgeColor = 'none';

% figure(3);
% plot(t,vse,'LineWidth',1.5);
% xlabel('time, ms');
% ylabel('Membrane potential [mV]');
% lgd=legend('0.25','0.5','0.75','1','1.25','1.5','1.75','2','2.25','2.5','2.75','3')
% title(lgd,'Neuron length, mm')
% grid on
% set(gca,'FontSize',15,'fontWeight','normal')

figure(3);
plot(t,vse,'LineWidth',1.5);
xlabel('Время, мс');
ylabel('Мембранный потенциал, мВ');
lgd=legend('0.25','0.5','0.75','1','1.25','1.5','1.75','2','2.25','2.5','2.75','3')
title(lgd,'Длина нейрона, мм')
grid on
set(gca,'FontSize',15,'fontWeight','normal')

figure(4);
s1=surf(l,t,vse);
xlabel('Length');
ylabel('Time');
zlabel('Voltage');
colorbar
s1.EdgeColor = 'none';

%  for j=1:100000
% 
%      view([j 25])
%      drawnow
%      pause(0.001)
%      
%  end

figure(5);
plot(r(l))



