clear all;

global vNa vK vL vH ...
    gNa gK gL gCa gSK gNap gKNa gH gA ...
    tauNbar tauHPbar tauE tauH Ca_ex RTF ...
    thetaM sigmaM thetaN sigmaN thetaS sigmaS ...
    thetaMP sigmaMP thetaHP sigmaHP thetaRF sigmaRF thetaRS sigmaRS ...
    kr prf taurs f ...
    thetaA sigmaA thetaE sigmaE ...
    C iapp eps kca ks pmax ec nh alphaNa rpump kp neq ...
    gT thetaaT sigmaaT thetabT sigmabT thetarT sigmarT taur0 taur1 thrT sgmrt phirT ...
    from to
  
% Nernest potentials in mV
vNa=50; 
vK=-90; % -85
vL=-70; % -65
vCa=80;
vH=-30;

% Conductances in nS
gNa=800;
gK=1700;
gL=2;
gCa=19;
gNap=1;
gSK=1;
gKNa=1;
gH=4;
gA=1;
gT=1.05;

% Time constants in mS
tauNbar=10;
tauHPbar=1000;
tauE=20;
tauH=1;
taurs=1500;
RTF=26.7;
Ca_ex=2.5;

% Function constants in mV
thetaM=-35;
thetaN=-30;
sigmaM=-5;
sigmaN=-5;
thetaS=-20;
sigmaS=-0.05;
thetaMP=-40;
sigmaMP=-6;
thetaHP=-48;
sigmaHP=6;
thetaA=-20;
sigmaA=-10;
thetaE=-60;
sigmaE=5;


thetaaT = -65;
sigmaaT = -7.8;
thetabT = 0.4;
sigmabT = -0.1;
thetarT = -67;
sigmarT = 2; 
taur0 = 40;
taur1 = 17.5;
thrT = 68;
sgmrt = 2.2;
phirT=0.2;

thetaRF=-105;
sigmaRF=5;
thetaRS=-105;
sigmaRS=25;

% Capacitance in pF
C=75;

% Other constants
kr=0.01;
prf=100;

f=0.1;
eps=0.0015;
kca=0.3;
ks=0.5;
pmax=0.37;
ec=38.7;
nh=3.5;
alphaNa=0.0001;
rpump=0.0008;
kp=15;
neq=8;

% Initial Conditons
v_i=-72.14; 
n_i=.0002;
ca_i=0.103;
h_i=0.01;
na_i=7.09;
rf_i=0.03;
rs_i=0.03;
e_i=0.05;
rT_i=0.01;
hp_i=0.1;

from=200;
to=1200;

j=1;

threshold = -30;

for i=40:20:200
    iapp=i;
   
    time=0:.01:1300;
    totalSize = size(time,2);
    
    init=[n_i;h_i;rf_i;rs_i;ca_i;v_i;na_i;e_i;rT_i;hp_i];
    
    [time,output]=ode113('DiffEquations',time,init);
    
    v = output(:,6);
    
    numOfSpikes=0;
    
        k=1;
        while(k<=totalSize)
            if(v(k)>=threshold)
                while(k<=totalSize && v(k)>=threshold)
                    k=k+1;
                end            
                numOfSpikes = numOfSpikes + 1;
            else
                k=k+1;
            end        
        end        
    
     totalNumOfSpikes(j) = numOfSpikes;
     
     frequency(j) = ((numOfSpikes/1000)*1000);
     x(j)=i;
     spikes(j)= numOfSpikes;
     j=j+1;  
              
end


figure();
plot(x,frequency, 'kv', 'Linewidth', 2);
ylabel('Average Firing Frequency (Hz)');
xlabel('Injected Current (I_{app})');
set(gca,'FontWeight','bold');
set(gca, 'Box', 'off' );
yl = get(gca, 'ylabel');
set(yl, 'FontWeight', 'bold') 
set(yl, 'FontSize', 14);
xl = get(gca, 'xlabel');
set(xl, 'FontWeight', 'bold') 
set(xl, 'FontSize', 14);
t = get(gca, 'title');
set(t, 'FontWeight', 'bold') 
set(t, 'FontSize', 14);
set(gca,'FontWeight','bold');
set(gca, 'Box', 'off' );

