%close all;
clear all;


%simulation time
T_MAX = 1000;
dt = 2.5e-3;

%List of time points
T = linspace(0, T_MAX, T_MAX/dt+1);



%create vectors of membrane terms over time

%voltage
v = zeros(length(T),1); 

%gating variables
mi = zeros(length(T),1);
hi = zeros(length(T),1); 
md = zeros(length(T),1);
mb = zeros(length(T),1);
mk = zeros(length(T),1);
ma = zeros(length(T),1); 
ha = zeros(length(T),1); 

%membrane currents
Ii = zeros(length(T),1); 
Id = zeros(length(T),1); 
Ib = zeros(length(T),1); 
Ik = zeros(length(T),1); 
Ia = zeros(length(T),1); 
Ic = zeros(length(T),1); 
Il = zeros(length(T),1); 



%initialize membrane variables.  These come from running model with
%I=-2+sin(2*pi/1000*T) (which induces burst-like firing) and no IC bursting current, for 50000 time units,
%and taking values at "final" timestep.  mb variable is one that takes a
%while to come to equilibrium

v(1) = -26.9349;

mi(1) = 0.3695;
hi(1) = 0.0226;
md(1) = 0.1174;
mb(1) = 0.6581;
mk(1) = 0.3415;
ma(1) = 0.9216;
ha(1) = 0.0033;

Ii(1) = -2.4375;
Id(1) = -4.7768;
Ib(1) = -1.9151;
Ik(1) = 7.4585;
Ia(1) = 1.0135;
Ic(1) = 0;
Il(1) = -1.3233;


%additional current that neuron receives
%Uncomment one of these options and explore the resulting dynamics!

%This should produce steady ("tonic") firing with no IC
I=-1*ones(length(T),1);

%This should produce bursts of firing with no IC
%I = -2+sin(2*pi/1000*T);

%This should produce no firing with no IC
%I=-3*ones(length(T),1);



vi = 45;
vd = 130;
vb = 85;
vk = -52;
va = -63;
vc = -73;
vl = 7;


%gating variables
%i
tau_mi = @(v) (11.5/(1+exp((v+16)/5.75)))+0.5;
mi_inf = @(v) (1 / (1+exp((v- -30.1)/-7)))^(3*2.7);
tau_hi = @(v) (85/(1+exp((v+24)/3.75)))+4.5;
hi_inf = @(v) (1 / (1+exp((v+19.5)/4.5)))^6; 

%d
tau_md = @(v) 7;
md_inf = @(v) (1 / (1+exp((v+3.5)/-7.25)))^1.7;

%b
tau_mb = @(v) (3000/(1+exp((v+40)/8.5)))+300;
mb_inf = @(v) 1+v/60;

%k
tau_mk = @(v) 32/(1+exp((v+18)/4.5))+15;
mk_inf = @(v) (1 / (1+exp((v- -27.6)/-14)))^(2*2.3);

%a
tau_ma = @(v) 18;
ma_inf = @(v) (1 / (1+exp((v+59)/-10.5)))^4;
tau_ha = @(v) 260;
ha_inf = @(v) (1 / (1+exp((v+73.9)/4)))^0.5;

%c, not included


%l, leak nothing to define here


%gate Ii conductance
gi0 = 29.7;
gi = gi0;


%for spike detection
thresh=0;
numspike=1;
spiketime=[0];

%membrane Cs is 19.8 nF
%do euler method for solving hh-type equations

for t = 1:length(T)-1
    
    %NO C
    v(t+1) = v(t) + (dt/(19.8e-3))*(-Il(t) -Ii(t) -Id(t) -Ib(t) -Ik(t) -Ia(t) +I(t) - Ic(t)) ;

    %gating variables
    mi(t+1) = mi(t) + dt*(tau_mi(v(t))^-1*(mi_inf(v(t)) - mi(t)));
    hi(t+1) = hi(t) + dt*(tau_hi(v(t))^-1*(hi_inf(v(t)) - hi(t)));   
    md(t+1) = md(t) + dt*(tau_md(v(t))^-1*(md_inf(v(t)) - md(t)));
    mb(t+1) = mb(t) + dt*(tau_mb(v(t))^-1*(mb_inf(v(t)) - mb(t)));
    mk(t+1) = mk(t) + dt*(tau_mk(v(t))^-1*(mk_inf(v(t)) - mk(t)));
    ma(t+1) = ma(t) + dt*(tau_ma(v(t))^-1*(ma_inf(v(t)) - ma(t)));
    ha(t+1) = ha(t) + dt*(tau_ha(v(t))^-1*(ha_inf(v(t)) - ha(t))); 
    
  

    %conductance gating
% Commented out.  Makes little qual. difference on test tims and introduces
% some oscillatory spike terms in Euler method with reasonable timesteps
%     if v(t+1) > 0
%         gi = 0;
%     else
%         gi = gi0;
%     end

    

    %currents
    Ii(t+1) = gi*(mi(t)^3*hi(t))*(v(t)-vi);
    Id(t+1) = .259*(md(t))*(v(t)-vd);
    Ib(t+1) = .026*(mb(t))*(v(t)-vb);
    Ik(t+1) = 2.55*(mk(t))^2*(v(t)-vk);
    Ia(t+1) = 11.7*(ma(t)^4*ha(t))*(v(t)-va); %inactivation above -30mv
    Il(t+1) = 0.039*(v(t)-vl);
    
          %condition for spike detection:  decreasing now, increased
          %before, over threshold
          if ((t>1) & (v(t+1)<v(t)) & (v(t)>v(t-1))) & v(t)>thresh & T(t)>(spiketime(numspike)+10)  %last condition ensures no two spikes within 10 ms
            numspike = numspike + 1;
            spiketime(numspike) = T(t+1);     
            
          end
    
    
    %display progress through simulation
    if mod(t,10000)==0
        disp(sprintf('fraction through simulation time = %g',t/length(T)))
    end

end



%%





%Plot voltage, current, gating variable traces vs time

figure
hold on;

%plot(T, mi, 'b');
%plot(T, hi, 'r');
%plot(T, md, 'c');

subplot(3,1,1);
plot(T, v, 'r');
ylabel('Voltage (mV)');
hold on
%add * at spike times at 0 voltage
plot(spiketime,0,'*')



subplot(3,1,2);
hold on;
plot(T, Ii, 'c');
plot(T, Id, 'g');
plot(T, Ib, 'b');
plot(T, Ik, 'k');
plot(T, Ia, 'm');
plot(T, Il, 'y');
ylabel('Current (nA)');



subplot(3,1,3);
hold on;
plot(T, mi, 'c');
plot(T, hi, 'c:');
plot(T, md, 'g');
plot(T, mb, 'b');
plot(T, mk, 'k');
plot(T, ma, 'm');
plot(T, ha, 'm:');

ylabel('Activation');
xlabel('Time (mSec)');

legend('mi','hi','md','mb','mk','ma','ha')