%% Fractional Leaky Integrate-and-Fire Model
% Model files for the paper
% Wondimu Teka,  Toma M.Marinov, Fidel Santamaria (2014) Neuronal Spike Timing Adaptation Described with a Fractional Leaky Integrate-and-Fire Model. 
% PLoS Comput Biol 10(3): e1003526. doi:10.1371/journal.pcbi.1003526

clear
 
rseed=2342;

Ncells=1; %number of cells,  do not change this value, this code is not for network, it is only for one cell. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %   parameters  values that can be varied. 

dt=0.1; %ms  time step. This value is used for  all the data. 
t=0:dt:1000; %ms    mostly is  up to  5000 ms (5sec)  or  10000 ms (10 sec) is used. 
                


 % alpha= ??; %  the fractional exponent (order)  it is  varied  below using a for loop

 
Rm=40; %Mohms
taum=20*ones([1 Ncells]); %ms membrane time constant
Cm=0.5; %nF  % this is not used since taum and Rm are used and note taum = Cm*Rm

refrac=8*ones([1 Ncells]); %ms absolute refarctory period 

v0= -70*ones([1 Ncells]);  % mV  initial value 
% initial voltage was varied: used values are v0= -70, v0= -75 v0= -58, v0= -45, 

vrest= -70;     % mV  the resting  potential which is used for reset too. 
vth= vrest + 20*ones([1 Ncells]); %  is  -50  mV,   threshold value of membrane voltage
 vpeak = vrest + 100*ones([1 Ncells]); % is  30 mV,   peak value of he voltage at spike


Iinjamplitude= 3;  % nA  injected current amplitude, see Iinj  for the whole injected current 

%  for the subthershold  use Iinjamplitude= 0.3; 


Namp=0; %noise amplitude, to add white Gaussian noise with this  standard devation Namp, which is noise amplitude

 %  This value should be changed only for figure 11, 
% Namp=1  %  For only figure 11,
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  %  used for most of the simulations
Iinj=Iinjamplitude*ones(length(t),Ncells)+ Namp*(randn(length(t),Ncells));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% for the sinusoidal simulation with  non constant freqency, use the following  ZAP current
%  f=0.1*(2./(1+exp(-t./1500))-1).^3;  % frequency  in cycle  per ms. since time is in ms; the freqency increases from 0  Hz to  100 Hz
%  Iinj=0.3*sin((2*pi).*f.*t)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %  for figure 10 use the following current instead of the above
% for the sinusoidal simulation with constant freqency, use the following
%f = 0.003  %  this is in  cycle per ms, it is the same is     3 Hz 
%Iinj= 2 +2*sin((0:dt:t(end)).*2*pi*f)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% names for parameters to pass them to the  main  function: runNetworkderivative

NetProp.Ncells=Ncells;
NetProp.Rm=Rm; %Mohms
NetProp.TauM=taum;
NetProp.Cm=Cm; %nF  % this is not used since taum and Rm are used.

NetProp.Refrac=refrac;
NetProp.vTh=vth;
NetProp.v0=v0;
NetProp.vrest=vrest;
NetProp.vpeak=vpeak;

NetProp.Noise=Namp;

NetProp.dt=dt;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%Main integrator
%  Integrating  output for different alpha values using  the fractional derivative

b=1;

 for alpha=1:-0.2:0.2
out(b)=runNetworkderivative(NetProp,Iinj,t,rseed,alpha); % main output
                % Note:This has the following output out.v=v  for voltage; out.sp=sp for spike ; out.Memo2=Memo2 for memery; out.t=t for time;

alphavalue(b)=alpha

  ISI{b}=diff(out(b).t(logical(out(b).sp(1:end-1))));  %  calculates  ISI in msec
  firingrate{b} = 1000./ISI{b};  %in  Hz, calculates firing rate
  totalspike(b) = sum(out(b).sp); % for toltal number of spikes for  total time
  b= b+1
 end
 
  cv=['kbgmr'];
  c=1;
 
  for b=1:length(out)  % for  alpha =1.0
  subplot(5,1,c)

      trange=1:1000/dt;
 
  plot(out(b).t(trange), out(b).v(trange), cv(b), 'Linewidth',1.5)  %


set(gca, 'FontSize', 14, 'FontName', 'Helvetica')

 set(gca, 'LooseInset', get(gca,'TightInset'))
 
 
ylabel('V (mV)', 'fontsize',14, 'FontName', 'Helvetica');
  legendinfo2= sprintf('\\alpha= %0.1f\n', alphavalue(b));
  c=c+1;
   axis([-1  1000  -75   35])
   
   legend(legendinfo2,'Location','SouthWest', 'fontsize',11, 'FontName', 'Helvetica')
   clear legendinfo2;
end

  xlabel('Time (ms)', 'fontsize',14, 'FontName', 'Helvetica');