%% this funcion is used in fractional neuron integration. It integrates the fractional derivative and  the  voltage v at each time t.

function out=runNetworkderivativeHHFractionalNa_h(NetProp,Iinj,t,alpha)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ncells=NetProp.Ncells;
dt=NetProp.dt;
Cm=NetProp.Cm;
v0=NetProp.v0;
vrest=NetProp.vrest;
gK=NetProp.gK;
gNa=NetProp.gNa;
gL=NetProp.gL;
EK=NetProp.EK;
ENa=NetProp.ENa;
EL=NetProp.EL;
m=NetProp.m;
h=NetProp.h;
n=NetProp.n;

Namp=NetProp.Noise;

%=================================
% variables in vecotr form 
v=vrest.*ones(length(t),Ncells);

mV=m*ones(length(t),Ncells);
hV=h*ones(length(t),Ncells);
nV=n*ones(length(t),Ncells);

VMemory=zeros(length(t)-1,Ncells);
Ngatememory=zeros(length(t)-1,Ncells);
Hgatememory=zeros(length(t)-1,Ncells);

I_HionicV=0*ones(length(t),Ncells);
INaV=0*ones(length(t),Ncells);
IKV=0*ones(length(t),Ncells);
ILV=0*ones(length(t),Ncells);


%============================================================
% The weight for the   memory trace of the fractional drivative for
% calculated here for the total time t for faster simulation
NN=length(t);
nn=1:NN-1;
WCoet=(NN+1-nn).^(1-alpha)-(NN-nn).^(1-alpha);

% Iionic=@(v,m,h,n,gbarL,gbarNa,gbarK,vrest,Ena,Ek)((gL*(v-EL)+...
%         gNa*m^3*h*(v-Ena)+...
%         gK*n^4*(v-Ek)));
   
       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
preT=1;

for a=preT
    [v(a+1,1),m,h,n]=HH_RKfull(v(a,1),m,h,n,gL,gNa,gK,EL,ENa,EK,Iinj(a),Cm,dt,v0);
   
    mV(a+1,:) = m;
    hV(a+1,:) = h;
    nV(a+1,:) = n;

     VMemory(a,:)=0;
    Ngatememory(a,:)=0;
    
       Hgatememory(a,:)=0;
end
%========================================================================================

kr = dt^alpha*gamma(2-alpha);     %  the kernel   from the fractional derivative that  weighted  the markovian term

for a=(preT(end)+1):length(t)-1   
    
   %Iionic =gL*(v(a,1)-EL)+gNa*m^3*h*(v(a,1)-ENa)+gK*n^4*(v(a,1)-EK);
           
    [v(a+1,1),m,n]=HH_RK(v(a,1),m,h,n,gL,gNa,gK,EL,ENa,EK,Iinj(a),Cm,dt,v0);
   
    
    
   
    
    mV(a+1,:) = m;
    nV(a+1,:) = n;
    %nV(a+1,:) = n;
    
    
       
    
    %%%%% The weight of the memory trace  
    WCoe=WCoet(end-a+2:end);  % The weight for the  voltage memory trace of the fractional drivative  at each  tiime t
    
    
  % voltage memory trace  
%     TeDi=v(2:a,:)-v(1:a-1,:); % Delta V (using all past values of V)  of  the  voltage memory trace of the fractional drivative  at each  tiime t
%     VoltMemory=WCoe*TeDi;  %  voltage memory trace 
%      VMemory(a,:)= VoltMemory;
%     fraccalcu=VoltMemory-v(a,1);  %  The fraction derivative
    
    %v(a+1,:) = (kr/Cm)*(Iinj(a)- Iionic) + v(a,1) - VoltMemory;
    
   %  v(a+1,:) = (dt/Cm)*(Iinj(a)- Iionic) + v(a,1) ;
     
     
    
    %%% Memory trace for  for gating variable n
    
%    DeltaN =nV(2:a,:)-nV(1:a-1,:);
%    NgateMemory=WCoe*DeltaN;
%    Ngatememory(a,:)= NgateMemory;
%       
%    alphan=(0.1-0.01*(v(a,1)-v0))./(exp(1-0.1*(v(a,1)-v0))-1);
%    betan=0.125.*exp(-(v(a,1)-v0)./80);
%   
%    
%     n = kr*(alphan*(1-nV(a,1))-betan*nV(a,1))+nV(a,1)-NgateMemory;
%    
%     nV(a+1,:) = n;
   


   DeltaH =hV(2:a,:)-hV(1:a-1,:);
   HgateMemory=WCoe*DeltaH;
   Hgatememory(a,:)= HgateMemory;
   
   alphah=0.07*exp(-(v(a,1)-v0)/20);
   betah=1./(exp(3-0.1*(v(a,1)-v0))+1);
      
  
  
   
    h = kr*(alphah*(1-hV(a,1))-betah*hV(a,1))+hV(a,1)-HgateMemory;
   
    hV(a+1,:) = h;



    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   

end
out.v=v;
%out. VMemory= VMemory;
out.t=t;
%out.fraccalcu=fraccalcu;
%out.I_Hionic=0;%I_Hionic;
%out.INa=INa;
%out.IK=IK;
%out.IL=IL;
out.mV=mV;
out.nV=nV;
out.hV=hV;

out.Ngatememory= Ngatememory;
out.Hgatememory= Hgatememory;
end


%function [vnew,mnew,hnew]=HH_RK(vold,mold,h1old,nold,dt,Iapp)
function [vnew,mnew,nnew]=HH_RK(vold,mold,h1old,nold,gL,gNa,gK,EL,ENa,EK,Iapp,Cm,dt,v0)


%v0=-65;
V=vold;
alpham=(2.5-0.1.*(V-v0))./(exp(2.5-0.1.*(V-v0))-1);
betam=4.*exp(-(V-v0)./18);
mf=@(Mm)(alpham*(1-Mm)-betam*Mm);


alphan=(0.1-0.01*(V-v0))./(exp(1-0.1*(V-v0))-1);
betan=0.125.*exp(-(V-v0)./80);
nf=@(Nna)(alphan*(1-Nna)-betan*Nna);

 


%%%%%%%%
k1=mf(mold);
k2=mf(mold+(dt/2)*k1);
k3=mf(mold+(dt/2)*k2);
k4=mf(mold+dt*k3);
mnew=mold+(dt/6)*(k1+2*k2+2*k3+k4);


k1=nf(nold);
k2=nf(nold+(dt/2)*k1);
k3=nf(nold+(dt/2)*k2);
k4=nf(nold+dt*k3);
nnew=nold+(dt/6)*(k1+2*k2+2*k3+k4);




% k1=hf(h1old);
% k2=hf(h1old+(dt/2)*k1);
% k3=hf(h1old+(dt/2)*k2);
% k4=hf(h1old+dt*k3);
% hnew=h1old+(dt/6)*(k1+2*k2+2*k3+k4);

%%%%%%%%%%%

%gK=36; gNa=120; gL=0.3;      % channel conductances: mS/cm2
%EK=-12 + v0(1,1); ENa=115 + v0(1,1); EL=10.6 + v0(1,1); % channel reversal potentials: mV
%Cm=1;

f=@(v,m,h,n,I)(-(1/Cm)*(gL*(v-EL)+gNa*m^3*h*(v-ENa)+gK*n^4*(v-EK))+I);
    
 
k1=f(vold,mold,h1old,nold,Iapp);
k2=f(vold+(dt/2)*k1,mold,h1old,nold,Iapp);
k3=f(vold+(dt/2)*k2,mold,h1old,nold,Iapp);
k4=f(vold+dt*k3,mold,h1old,nold,Iapp);
vnew=vold+(dt/6)*(k1+2*k2+2*k3+k4);

end





function [vnew,mnew,hnew,nnew]=HH_RKfull(vold,mold,h1old,nold,gbarL,gbarNa,gbarK,EL,Ena,Ek,Im,Cm,dt,v0)


V=vold;
alpham=(2.5-0.1.*(V-v0))./(exp(2.5-0.1.*(V-v0))-1);
betam=4.*exp(-(V-v0)./18);
alphah=0.07*exp(-(V-v0)/20);
betah=1./(exp(3-0.1*(V-v0))+1);
alphan=(0.1-0.01*(V-v0))./(exp(1-0.1*(V-v0))-1);
betan=0.125.*exp(-(V-v0)./80);
mf=@(Mm)(alpham*(1-Mm)-betam*Mm);
hf=@(Hh)(alphah*(1-Hh)-betah*Hh);
nf=@(Nna)(alphan*(1-Nna)-betan*Nna);



f=@(v,m,h,n,gbarL,gbarNa,gbarK,vrest,Ena,Ek,I,Cm)(-(1/Cm)*(gbarL*(v-EL)+...
        gbarNa*m^3*h*(v-Ena)+...
        gbarK*n^4*(v-Ek))+...
        I);
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k1=f(vold,mold,h1old,nold,gbarL,gbarNa,gbarK,EL,Ena,Ek,Im,Cm);
k2=f(vold+(dt/2)*k1,mold,h1old,nold,gbarL,gbarNa,gbarK,EL,Ena,Ek,Im,Cm);
k3=f(vold+(dt/2)*k2,mold,h1old,nold,gbarL,gbarNa,gbarK,EL,Ena,Ek,Im,Cm);
k4=f(vold+dt*k3,mold,h1old,nold,gbarL,gbarNa,gbarK,EL,Ena,Ek,Im,Cm);
vnew=vold+(dt/6)*(k1+2*k2+2*k3+k4);


k1=mf(mold);
k2=mf(mold+(dt/2)*k1);
k3=mf(mold+(dt/2)*k2);
k4=mf(mold+dt*k3);
mnew=mold+(dt/6)*(k1+2*k2+2*k3+k4);

k1=hf(h1old);
k2=hf(h1old+(dt/2)*k1);
k3=hf(h1old+(dt/2)*k2);
k4=hf(h1old+dt*k3);
hnew=h1old+(dt/6)*(k1+2*k2+2*k3+k4);


k1=nf(nold);
k2=nf(nold+(dt/2)*k1);
k3=nf(nold+(dt/2)*k2);
k4=nf(nold+dt*k3);
nnew=nold+(dt/6)*(k1+2*k2+2*k3+k4);


end