% HGR Apr05-13
% HGR Mar23-13

close all;
clearvars;

MDL = 1;
         % 1: hei-vav-pei-fei
         % 2: samech-nun-iod-c        
         % 3: hei-vav-mem-vav-kuf-lamed-iod-nun-iod-c
KSE = 1;
         % 1: Single cell
         % 2: Two-cell network

         
% Biophysical Parameters

C = 20;
El = -60;
Eca = 120;
Ek = -84;
Gl = 2;
Gk = 6;
V1 = -1.2;
V2 = 18;
if MDL == 1
    Gca = 4.4;
    Phi = 0.04;
    V3 = 2;
    V4 = 30;
    Iapp = 80;
elseif MDL == 2
    Gca = 4;
    Phi = 0.067;
    V3 = 12;
    V4 = 17.4;
    Iapp = 40;
elseif MDL == 3
    Gca = 4;
    Phi = 0.23;
    V3 = 12;
    V4 = 17.4;
    Iapp = 30;
end

Phi = 0.01;

Gex = 0.4;
Gin = 0;
MVN = 0;
Iapp = 80;
t1appl = 0;
t2appl = 150;
tconn = 1000;
tconnoff = 4000;


% Activation/inactivation curves and time constant

minf=@(v) 0.5 * (1 + tanh((v - V1) / V2));
winf=@(v) 0.5 * (1 + tanh((v - V3) / V4));
tauw=@(v) Phi * (cosh((v - V3) / (2 * V4)));

% Sigmoid function

%vsth = 0;
%vsslp = 4;
vsth = 2;
vsslp = 5;
H=@(v) 0.5*(1 + tanh(v-vsth/vsslp)); 


Ain = 0;
f  = 50;


if KSE == 1
    % Nullclines

    vv=-100:0.1:100;
    vnlc = (Iapp - Gl*(vv-El)-Gca*minf(vv).*(vv-Eca))./(Gk*(vv-Ek));
    

    % Numerical Solution


    T = 2000;
    dt = 0.1;
    t = 0:dt:T;

    V = zeros(1,length(t));
    w = zeros(1,length(t));

    V(1) = -60;
    w(1) = 0.4;

    Iin = Ain*sin(2*pi*f*t/1000);
    
    for j=1:length(t)-1
        k1v = Iapp-Gl*(V(j)-El)-Gk*w(j)*(V(j)-Ek)-Gca*minf(V(j))*(V(j)-Eca);
        k1v = k1v + Iin(j);
        k1v = k1v / C;
        k1w = (winf(V(j)) - w(j))*tauw(V(j));
        av = V(j) + k1v*dt;
        aw = w(j) + k1w*dt;
        k2v = Iapp-Gl*(av-El)-Gk*aw*(av-Ek)-Gca*minf(av)*(av-Eca);
        k2v = k2v + Iin(j+1);
        k2v = k2v / C;
        k2w = (winf(av) - aw)*tauw(av);
        V(j+1) = V(j) + (k1v + k2v)*dt/2;
        w(j+1) = w(j) + (k1w + k2w)*dt/2;
    end

    Tmin=T-1000;
    Tmax=T;
    [spkt spkcnt]= spiketimes(V,t,Tmin,Tmax);
    ISI=diff(spkt(1:spkcnt));

    figure
    plot(t,V,'linewidth',2);
    axis([0 Tmax -80 60]);
    set(gca,'fontsize',20);
    xlabel('t');
    ylabel('V');

    figure
    plot(vv,vnlc,'r','linewidth',2);
    hold on
    plot(vv,winf(vv),'g','linewidth',2);
    plot(V,w,'linewidth',2);
    plot(V(1),w(1),'ob','linewidth',3);
    axis([-80 100 -0.1 1]);
    set(gca,'fontsize',20);
    xlabel('V');
    ylabel('w');
    legend('v-nullcline','w-nullcline','trajectory','IC')

    MVN = 0;
    if MVN == 1
        figure
        k = 50;
        for j=1:floor(length(t)/k)-1
            vnlcfrzd = vnlc + Iin(k*j)./(Gk*(vv-Ek));
            plot(vv,vnlc,'--r','linewidth',1);
            hold on
            plot(vv,winf(vv),'g','linewidth',2);
            plot(vv,vnlcfrzd,'-r','linewidth',2);
            plot(V(k*j:k*j),w(k*j:k*j),'o','linewidth',2);
            plot(V(1),w(1),'ob','linewidth',3);
            axis([-80 100 -0.1 1]);
            set(gca,'fontsize',20);
            xlabel('V');
            ylabel('w');
            legend('v-nullcline','w-nullcline','trajectory','IC')
            pause(0.0001)
            hold off
        end
    end
elseif KSE == 2
    % Numerical Solution

    taurse_e = 1;     % AMPA
    taudec_e = 5;     % AMPA
    taurse_i = 0.2;     % GABA_A     
    taudec_i = 10;     % GABA_A
    Eex = 0;
    Ein = -30;

    %Iapp = 80;
    %Gex = 0;
    %Gin = 0.5;
    Ggap = 0;

   
    T = 5000;
    dt = 0.1;
    t = 0:dt:T;
%     t1appl = 0;
%     t2appl = 10;
%    tconn = 1000;
  
    
    V1 = zeros(1,length(t));
    w1 = zeros(1,length(t));
    Se1 = zeros(1,length(t));
    Si1 = zeros(1,length(t));
    V2 = zeros(1,length(t));
    w2 = zeros(1,length(t));
    Se2 = zeros(1,length(t));
    Si2 = zeros(1,length(t));

    IC = 3;
    
    if IC == 1
        V1(1) = -60;
        w1(1) = 0.6;
        V2(1) = -59.5;
        w2(1) = 0.6; 
    elseif IC == 2
        V1(1) = -50;
        w1(1) = 0.4;
        V2(1) = 40;
        w2(1) = 0.4; 
    elseif IC == 3
        V1(1) = -60;
        w1(1) = 0;
        V2(1) = V1(1);
        w2(1) = w1(1);        
    end
    
    Se1(1) = 0;
    Si1(1) = 0;
    Se2(1) = 0;
    Si2(1) = 0;

   
   
   
    for j=1:length(t)-1
        if t(j) < tconn
            Gex21 = 0;
            Gex12 = 0;
            Gin12 = 0;
            Gin21 = 0;
            Gc = 0;
        else
            Gex21 = Gex;
            Gex12 = Gex;
            Gin12 = Gin;
            Gin21 = Gin;
            Gc = Ggap;
        end
        if t(j)>tconnoff
            Iapp = 60;
        end
        
        k1v1 = Iapp-Gl*(V1(j)-El)-Gk*w1(j)*(V1(j)-Ek)-Gca*minf(V1(j))*(V1(j)-Eca);       
        k1v1 = k1v1 + Ain*sin(2*pi*f*t(j)/1000); 
        k1v1 = k1v1 - Gex21*Se2(j)*(V1(j)-Eex);
        k1v1 = k1v1 - Gin21*Si2(j)*(V1(j)-Ein);
        k1v1 = k1v1 + Gc*(V2(j)-V1(j));
        k1v1 = k1v1 / C;
        k1w1 = (winf(V1(j)) - w1(j))*tauw(V1(j));       
        k1se1 = H(V1(j))*(1-Se1(j))/taurse_e-Se1(j)/taudec_e;       
        k1si1 = H(V1(j))*(1-Si1(j))/taurse_i-Si1(j)/taudec_i;
        k1v2 = Iapp-Gl*(V2(j)-El)-Gk*w2(j)*(V2(j)-Ek)-Gca*minf(V2(j))*(V2(j)-Eca);
        k1v2 = k1v2 + Ain*sin(2*pi*f*t(j)/1000);
        k1v2 = k1v2 - Gex12*Se1(j)*(V2(j)-Eex);
        k1v2 = k1v2 - Gin12*Si1(j)*(V2(j)-Ein);
        k1v2 = k1v2 + Gc*(V1(j)-V2(j));
        k1v2 = k1v2 / C;
        k1w2 = (winf(V2(j)) - w2(j))*tauw(V2(j));
        k1se2 = H(V2(j))*(1-Se2(j))/taurse_e-Se2(j)/taudec_e;
        k1si2 = H(V2(j))*(1-Si2(j))/taurse_i-Si2(j)/taudec_i;
        av1 = V1(j) + k1v1*dt;
        aw1 = w1(j) + k1w1*dt;
        ase1 = Se1(j) + k1se1*dt;
        asi1 = Si1(j) + k1si1*dt;
        av2 = V2(j) + k1v2*dt;
        aw2 = w2(j) + k1w2*dt;
        ase2 = Se2(j) + k1se2*dt;
        asi2 = Si2(j) + k1si2*dt;
        k2v1 = Iapp-Gl*(av1-El)-Gk*aw1*(av1-Ek)-Gca*minf(av1)*(av1-Eca);
        k2v1 = k2v1 + Ain*sin(2*pi*f*t(j+1)/1000);
        k2v1 = k2v1 - Gex21*ase2*(av1-Eex);
        k2v1 = k2v1 - Gin21*asi2*(av1-Ein);
        k2v1 = k2v1 + Gc*(av2-av1);
        k2v1 = k2v1 / C;
        k2w1 = (winf(av1) - aw1)*tauw(av1);
        k2se1 = H(av1)*(1-ase1)/taurse_e-ase1/taudec_e;
        k2si1 = H(av1)*(1-asi1)/taurse_i-asi1/taudec_i;
        k2v2 = Iapp-Gl*(av2-El)-Gk*aw2*(av2-Ek)-Gca*minf(av2)*(av2-Eca);
        k2v2 = k2v2 + Ain*sin(2*pi*f*t(j+1)/1000);
        k2v2 = k2v2 - Gex12*ase1*(av2-Eex);
        k2v2 = k2v2 - Gin12*asi1*(av2-Ein);
        k2v2 = k2v2 + Gc*(av1-av2);
        k2v2 = k2v2 / C;
        k2w2 = (winf(av2) - aw2)*tauw(av2);
        k2se2 = H(av2)*(1-ase2)/taurse_e-ase2/taudec_e;
        k2si2 = H(av2)*(1-asi2)/taurse_i-asi2/taudec_i;
        V1(j+1) = V1(j) + (k1v1 + k2v1)*dt/2;
        w1(j+1) = w1(j) + (k1w1 + k2w1)*dt/2;
        V2(j+1) = V2(j) + (k1v2 + k2v2)*dt/2;
        w2(j+1) = w2(j) + (k1w2 + k2w2)*dt/2;
        Se1(j+1) = Se1(j) + (k1se1 + k2se1)*dt/2;
        Se2(j+1) = Se2(j) + (k1se2 + k2se2)*dt/2;
        Si1(j+1) = Si1(j) + (k1si1 + k2si1)*dt/2;
        Si2(j+1) = Si2(j) + (k1si2 + k2si2)*dt/2;
        if t(j) < t1appl
            V1(j+1) = V1(j);
            w1(j+1) = w1(j);
        end
        if t(j) < t2appl
            V2(j+1) = V2(j);
            w2(j+1) = w2(j);
        end
        
    end

    Tmin=500;
    Tmax=T;
    [spkt1 spkcnt1]= spiketimes(V1,t,Tmin,Tmax);
    ISI1=diff(spkt1(1:spkcnt1));
    [spkt2 spkcnt2]= spiketimes(V2,t,Tmin,Tmax);
    ISI2=diff(spkt2(1:spkcnt2));
    if spkcnt1 > 2 && spkcnt2 > 2
        Phase = abs(spkt2(spkcnt2-2)-spkt1(spkcnt1-2))/(spkt1(spkcnt1-1)-spkt1(spkcnt1-2));
        if Phase > 0.5
            Phase = 1-Phase;
        end

        for j=1:spkcnt1-2
            Phse(j) = abs(spkt2(j+1)-spkt1(j+1))/(spkt1(j+1)-spkt1(j));
            if Phse(j) > 0.5
                Phse(j) = 1 - Phse(j);
            end
            if spkt1(j+1)<tconn && spkt2(j+1)<tconn
                jphseconn = j;
            end
        end
    end
    
    
    
    figure
    plot(-1000,-1000,'b','linewidth',2);
    hold on
    plot(-1000,-1000,'r','linewidth',2);
    plot([tconn tconn],[-100 100],'--','Color',[.7 .7 .7],'linewidth',2);
    plot([tconnoff tconnoff],[-100 100],'--','Color',[.7 .7 .7],'linewidth',2);
    plot(t,V1,'b','linewidth',2);
    hold on
    plot(t,V2,'r','linewidth',2);
    axis([0 Tmax -80 80]);
    set(gca,'fontsize',24);
    xlabel('t  [ms]');
    ylabel('V  [mV]');
    %legend('V_1','V_2');
    
    if spkcnt1 > 2 && spkcnt2 > 2
        figure
        plot(-100 -100,'-o','linewidth',2);
        hold on
        plot([0 100],[0 0],'--','Color',[.5 .5 .5],'linewidth',1);
        plot([0 100],[0.5 0.5],'--','Color',[.5 .5 .5],'linewidth',1);
        plot([jphseconn jphseconn],[-1 1],'--','Color',[.7 .7 .7],'linewidth',1.5);
        plot(Phse,'-o','linewidth',2)
        axis([0 length(Phse) -0.05 0.55]);
        set(gca,'fontsize',24);
        xlabel('cycle number');
        ylabel('\Phi');
    end
    
    %MVN = 0;
    
    if MVN == 1
        % Nullclines

        vv=-100:0.1:100;
        vnlc = (Iapp - Gl*(vv-El)-Gca*minf(vv).*(vv-Eca))./(Gk*(vv-Ek));
        
        figure
        k = 10;
        for j=floor((tconn-850)/dt/k):floor(length(t)/k)-1
            vnlc1 = vnlc - (Gex21*Se2(k*j)*(V1(k*j)-Eex)+Gin21*Si2(k*j)*(V1(k*j)-Ein))./(Gk*(vv-Ek));
            vnlc2 = vnlc - (Gex12*Se1(k*j)*(V2(k*j)-Eex)+Gin12*Si1(k*j)*(V2(k*j)-Ein))./(Gk*(vv-Ek));
            plot(vv,vnlc,'--','Color',[.7 .7 .7],'linewidth',1.5);
            hold on
            plot(vv,winf(vv),'g','linewidth',2);
            plot(vv,vnlc1,'b','linewidth',2);
            plot(vv,vnlc2,'r','linewidth',2);
            plot(V1(k*j:k*j),w1(k*j:k*j),'ob','linewidth',2);
            plot(V2(k*j:k*j),w2(k*j:k*j),'or','linewidth',2);
            axis([-60 60 0 0.7]);
            set(gca,'fontsize',20);
            xlabel('V');
            ylabel('w');
            pause(0.0001)
            hold off
        end
    end
    
    
end