%% This is a funcion for the Fractional Leaky IntegrateandFire Model. It integrates the fractional derivative and the voltage v at each time t.
function out=runNetworkderivative(NetProp,Iinj,t,rseed,alpha)
rand('seed',rseed)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ncells=NetProp.Ncells;
taum=NetProp.TauM;
refrac=NetProp.Refrac/NetProp.dt;
vth=NetProp.vTh;
vpeak=NetProp.vpeak;
v0=NetProp.v0;
vrest=NetProp.vrest;
Namp=NetProp.Noise;
rm=NetProp.Rm;
dt=NetProp.dt;
v=vrest.*ones(length(t),Ncells);
v(1,:) = v0(1,:);
vrest=vrest.*ones(length(t),Ncells);
sp=zeros(length(t),Ncells);
isstillrefrac=zeros(1,Ncells);
rfcounter=zeros(1,Ncells);
Inhib = zeros(length(t),Ncells);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The weight for the voltage memory trace of the fractional drivative for
% calculated here for the total time t for faster simulation
NN=length(t);
nn=1:NN1;
WCoet=(NN+1nn).^(1alpha)(NNnn).^(1alpha);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for a=1:length(t)1
Ihere=Iinj(a,:);
if a<2 %for the first step, it is a normal neuron
vdummy =dt*(vrest(a,:)v(a,:)+rm*Ihere)./taum + v(a,:);
else
%%%%% Fractional derivative starts here
d2dM=v(1:a,:); % to call all past values of voltage for fractioanl integration
WCoe=WCoet(enda+2:end); % The weight for the voltage memory trace of the fractional drivative at each tiime t
TeDi=d2dM(2:a,:)d2dM(1:a1,:); % Delta V (using all past values of V) of the voltage memory trace of the fractional drivative at each tiime t
fraccalcu=WCoe*TeDid2dM(a,:); % The fraction derivative
kr = dt^alpha*gamma(2alpha); % the kernel from the fractional derivative and weighted the markovian term
%%%%% Fractional derivative ends here
% used for the fractioanl neuron, value of the voltage if it is on subthreshold.
vdummy = kr*(vrest(a,:)v(a,:)+rm*Ihere)./taum  fraccalcu;
Memo2(a,:)=fraccalcu+v(a,:); % to save memory over time
if Memo2(a,:) > 0 && Memo2(a1,:) < 0
Inhib(a1,:) = Memo2(a1,:); % to save the value of mmemory which has a positive feedback
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
test=((vdummy>vth)+(isstillrefrac==1))>0 ;%v larger than vth or inside refrac period
sp(a+1,:)=test.*(rfcounter==0); % a spike if outside refrac period and test is possitive
v(a+1,:)=vrest(a,:).*(isstillrefrac==1) + (~test).*vdummy + sp(a+1,:)*vpeak; % save V value
isstillrefrac=sp(a+1,:)+test.*isstillrefrac.*(rfcounter<refrac); %initiallize refrac flag
rfcounter=rfcounter+test.*(isstillrefrac==1); %you enter or continue rfcounter
rfcounter=(~(test.*(rfcounter==refrac))).*rfcounter; %exit rfcounter
isstillrefrac=isstillrefrac.*(rfcounter>0); %resetting the refrac flag
end
out.v=v;
out.sp=sp;
out.Memo2=Memo2;
out.t=t;
out.Inhib = Inhib;
