%%%% Three dimensional bifurcation code of Leech's heart interneuron model (See fig.1 in the manuscript).
%%%% Here I refers to V^{shift}_{K2} and u refers to the membrane voltage(v) of the manuscript.
%%To run this code put the peakdet.m and this file in the same folder.
%%% Title- Emergence of bursting in a network of memory dependent excitable and spiking Leech-Heart neurons
%%%  Authors-Sanjeev Kumar Sharma; Argha Monda; Arnab Mondal; Ranjit Kumar Upadhyay and Chittaranjan Hens


for alpha=0.7:0.01:1
    alpha
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 dt = 0.001; tspan = 0:dt:5;
 %%%%%%%% Initial conditions %%%%%%%%%%%%%%%%%%%%%%%%%
 u=rand/10;
 V=rand/10;
 w=rand/10;
%   u=(-0.0271517+rand); 
%   V=(0.0442+rand);
%   w=(0.0513+rand);
%   u=-0.0288274+rand;
%   V=0.0965+rand;
%   w=0.3067+rand;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
uu=zeros(length(tspan),1); VV=zeros(length(tspan),1); ww=zeros(length(tspan),1);
uu(1,1)=u;
VV(1,1)=V;
ww(1,1)=w;

T1=tspan(end)/10;
T2=T1+60;
NN=length(tspan);
nn=1:NN-1;

WCoet=(NN+1-nn).^(1-alpha)-(NN-nn).^(1-alpha);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for I=-0.04:0.001:0.04
     I
for j=1:length(tspan)-1    
%     if (tspan(j)>T1) 
%      I=42;
%      %I=-99;
%     else 
%         I=0;
%     end;
    if j<2
  u=u+dt*(-2*(30*(w.^2)*(u+0.07)+8*(u+0.046)+200*(1/(1+exp(-150*(0.0305+u))).^3)*V*(u-0.045)));
   V=V+dt*(24.69*(1/(1+exp(500*(0.0333+u)))-V));
  w=w+dt*(4*(1/(1+exp(-83*(0.018+I+u)))-w));
  
    else
  %%%%% Fractional derivative starts  here
       WCoe=WCoet(end-j+2:end);  % The weight   of the fractional drivative  at each  tiime t 
       kr = dt^alpha*gamma(2-alpha);     %  the kernel   from the fractional derivative and  weighted  the markovian term
        %%%%% Fractional derivative for u starts  here 
        
        d2dM=uu(1:j,:); % to call all past values of voltage u  for fractioanl integration
       
        TeDi=d2dM(2:j,:)-d2dM(1:j-1,:); % Delta uu (using all past values of u)  of  the  voltage memory trace of the fractional drivative  at each  tiime t 
        fraccalcu=WCoe*TeDi-d2dM(j,:);  %  The fraction derivative 
         %%%%% Fractional derivative ends   here 
        
        %%%%% Fractional derivative for V starts  here 
        d2dMV=VV(1:j,1); % to call all past values  VV  for fractioanl integration
        TeDiV=d2dMV(2:j,1)-d2dMV(1:j-1,1); % Delta V (using all past values of VV) 
        fraccalcuV=WCoe*TeDiV-d2dMV(j,1);  %  The fraction derivative 
        %%%%% Fractional derivative ends   here 
             %%%% ==w==  Fractional derivative for w starts  here 
        d2dMw=ww(1:j,1); % to call all past values ww  for fractioanl integration
        TeDiw=d2dMw(2:j,1)-d2dMw(1:j-1,1); % Delta w (using all past values of ww) 
        fraccalcuw=WCoe*TeDiw-d2dMw(j,1);  %  The fraction derivative 
        %%%%% Fractional derivative ends   here 
       
         u =kr*(-2*(30*(w.^2)*(u+0.07)+8*(u+0.046)+200*(1/(1+exp(-150*(0.0305+u))).^3)*V*(u-0.045)))- fraccalcu; 
         V =kr*(24.69*(1/(1+exp(500*(0.0333+u)))-V))-fraccalcuV;
         w=kr*(4*(1/(1+exp(-83*(0.018+I+u)))-w))-fraccalcuw;
        
%         Memo2(j,:)=fraccalcu+uu(j,:); % to save memory over time
%         Memo2V(j,:)=fraccalcuV+VV(j,:); % to save memory over time
%         Memo2w(j,:)=fraccalcuw+ww(j,:);
    end
  
    uu(j+1,1)=u;
     VV(j+1,1)=V;
     ww(j+1,1)=w;
end
y(:,1)=uu(:,1);
y(:,2)=VV(:,1);
y(:,3)=ww(:,1);
%%%%%%%%%%%%%%%%%%%%%
a=size(y(:,1))
b=round(3*a/4)
x=y(b:end,1)
[maxtab,mintab]=peakdet(x,0.001)
%%%%%%%%%%%%%%%%%%%%%%%%%% Plot command for different slices %%%%%%%%%%%%% 
if(alpha==0.82)
    
if((size(maxtab)~=[0 0]) & (size(mintab)~=[0 0]))
    set(gca,'FontName','Times New Roman','FontSize',24)
plot3(alpha,I,maxtab(:,2),'g .','markersize',2)
   xlabel(' \it \bf \alpha ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
ylabel(' \it \bf V^{shift}_{K2} ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
   zlabel(' \it \bf v ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
hold on
plot3(alpha,I,mintab(:,2) ,'g .','markersize',2)
 xlabel(' \it \bf \alpha ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
ylabel(' \it \bf V^{shift}_{K2} ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
   zlabel(' \it \bf v ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
else
     y1=y(size(y,1),:);
     y2=y1(:,1);
plot3(alpha,I,y2, 'g .','markersize',2)
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
if(alpha==0.87)
     
if((size(maxtab)~=[0 0]) & (size(mintab)~=[0 0]))
    set(gca,'FontName','Times New Roman','FontSize',24)
plot3(alpha,I,maxtab(:,2),'r .','markersize',2)
   xlabel(' \it \bf \alpha ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
ylabel(' \it \bf V^{shift}_{K2} ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
   zlabel(' \it \bf v ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
hold on
plot3(alpha,I,mintab(:,2) ,'r .','markersize',2)
 xlabel(' \it \bf \alpha ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
ylabel(' \it \bf V^{shift}_{K2} ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
   zlabel(' \it \bf v ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
else
     y1=y(size(y,1),:);
     y2=y1(:,1);
plot3(alpha,I,y2, 'r .','markersize',2)
end
end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(alpha==0.94)
    
if((size(maxtab)~=[0 0]) & (size(mintab)~=[0 0]))
    set(gca,'FontName','Times New Roman','FontSize',24)
plot3(alpha,I,maxtab(:,2),'b .','markersize',2)
   xlabel(' \it \bf \alpha ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
ylabel(' \it \bf V^{shift}_{K2} ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
   zlabel(' \it \bf v ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
hold on
plot3(alpha,I,mintab(:,2) ,'b .','markersize',2)
 xlabel(' \it \bf \alpha ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
ylabel(' \it \bf V^{shift}_{K2} ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
   zlabel(' \it \bf v ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
else
     y1=y(size(y,1),:);
     y2=y1(:,1);
plot3(alpha,I,y2, 'b .','markersize',2)
end

   end
%   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
   if(alpha==1)
    
if((size(maxtab)~=[0 0]) & (size(mintab)~=[0 0]))
    set(gca,'FontName','Times New Roman','FontSize',24)
plot3(alpha,I,maxtab(:,2),'k .','markersize',2)
   xlabel(' \it \bf \alpha ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
ylabel(' \it \bf V^{shift}_{K2} ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
   zlabel(' \it \bf v ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
hold on
plot3(alpha,I,mintab(:,2) ,'k .','markersize',2)
 xlabel(' \it \bf \alpha ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
ylabel(' \it \bf V^{shift}_{K2} ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
   zlabel(' \it \bf v ', 'FontName', 'Times New Roman', ...
       'FontSize',30,'Color','k', 'Interpreter', 'tex')
else
     y1=y(size(y,1),:);
     y2=y1(:,1);
plot3(alpha,I,y2, 'k .','markersize',2)
end

  end
 end
end

