Fractional leaky integrate-and-fire model (Teka et al. 2014)

 Download zip file 
Help downloading and running models
We developed the Fractional Leaky Integrate-and-Fire model that can produce downward and upward spike time adaptions observed on pyramidal cells.The adaptation emerges from the fractional exponent of the voltage dynamics.
1 . Teka W, Marinov TM, Santamaria F (2014) Neuronal spike timing adaptation described with a fractional leaky integrate-and-fire model. PLoS Comput Biol 10:e1003526 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Abstract integrate-and-fire fractional leaky neuron;
Channel(s): I Na, leak;
Gap Junctions:
Simulation Environment: MATLAB;
Model Concept(s): Spike Frequency Adaptation;
Search NeuronDB for information about:  I Na, leak;
%% 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


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.Rm=Rm; %Mohms
NetProp.Cm=Cm; %nF  % this is not used since taum and Rm are used.




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


 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;


  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
  for b=1:length(out)  % for  alpha =1.0

  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));
   axis([-1  1000  -75   35])
   legend(legendinfo2,'Location','SouthWest', 'fontsize',11, 'FontName', 'Helvetica')
   clear legendinfo2;

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


Loading data, please wait...