% 2020-01-05

% Nap moodel (1D)

clearvars;
close all;

% Functions

Vphlf = -38;
Vpslp = 6.5;
pinf=@(V)  1.0./(1+exp(-(V-Vphlf)/Vpslp));


% Biophysical parameters

C = 1;
El = -65;
Ena = 55;
Gl = 0.25;
Gp = 0.16;
D = 0;
Iapp = 0;

F=@(V) Iapp-Gl*(V-El)-Gp*pinf(V).*(V-Ena);



% Computation of the fixed-point (assumes at most 3 fixed-points, brute force)

vv=-100:0.001:100;

fpcnt  = 0;
if Gp == 0 && Gl > 0 
    Vfp = El+Iapp/Gl;
else
    if F(vv(1)) > 0
        for j=1:length(vv)
            if F(vv(j))<0 && fpcnt == 0
                Vfp(1) = (vv(j)+vv(j-1))/2;
                fpcnt = 1;
            elseif F(vv(j))>0 && fpcnt == 1
                Vfp(2) = (vv(j)+vv(j-1))/2;
                fpcnt = 2;
            elseif F(vv(j))<0 && fpcnt == 2
                Vfp(3) = (vv(j)+vv(j-1))/2;
                fpcnt = 3;
            end
        end
    end
end

% Computation of the eigenvalues (Note: they can be computed analytically)

dv = 0.001;
Eig = zeros(1,length(Vfp));
for k=1:length(Vfp)
    Eig(k) = (F(Vfp(k)+dv)-F(Vfp(k)-dv))/(2*dv);
end

% Phase-line diagram


figure(1)
hFig = figure(1);
set(hFig, 'Position', [40 400 1000 500]); 
subplot(1,2,1)
hold on
plot(vv,F(vv),'r','linewidth',2);
plot([-100 100],[0 0],'--','Color',[.7 .7 .7],'linewidth',1);
axis([-80 80 -10 30]);
set(gca,'fontsize',20);
xlabel('V  [mV]');
ylabel('dV / dt');
for k=1:length(Vfp)
    if Eig(k)<0
         plot(Vfp(k),0,'ok','linewidth',2);
         plot(Vfp(k),0,'xk','linewidth',2);
    elseif Eig(k)>0       
        plot(Vfp(k),0,'ok','linewidth',2);        
     end
end

% Computation of the numerical solutions


Tmax = 100;
dt = 0.1;
t = 0:dt:Tmax;
eta = randn(1,length(t));


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

if length(Vfp) == 1

    V(1) = -60;

    for j=1:length(t)-1  
        k1v = (Iapp-Gl*(V(j)-El)-Gp*pinf(V(j))*(V(j)-Ena))/C;                 
        av = V(j)+k1v*dt;
        av = av+sqrt(2*D*dt)*eta(j);    
        k2v = (Iapp-Gl*(av-El)-Gp*pinf(av)*(av-Ena))/C;     
        V(j+1) = V(j)+(k1v+k2v)*dt/2;
        V(j+1) = V(j+1)+sqrt(2*D*dt)*eta(j);
    end


    figure(1)
    subplot(1,2,2)
    hold on
    plot(t,V,'-b','linewidth',2);
    axis([0 Tmax -80 40]);
    set(gca,'fontsize',24);
    xlabel('t  [ms]');
    ylabel('V  [mV]');
    
    
    
elseif length(Vfp)== 3
    
    V(1) = Vfp(2)-0.1;

    for j=1:length(t)-1  
        k1v = (Iapp-Gl*(V(j)-El)-Gp*pinf(V(j))*(V(j)-Ena))/C;                 
        av = V(j)+k1v*dt;
        av = av+sqrt(2*D*dt)*eta(j);    
        k2v = (Iapp-Gl*(av-El)-Gp*pinf(av)*(av-Ena))/C;     
        V(j+1) = V(j)+(k1v+k2v)*dt/2;
        V(j+1) = V(j+1)+sqrt(2*D*dt)*eta(j);
    end

    V1 = V;

    figure(1)
    subplot(1,2,2)
    hold on
    plot(t,V,'-b','linewidth',2);
    axis([0 Tmax -80 40]);
    set(gca,'fontsize',24);
    xlabel('t  [ms]');
    ylabel('V  [mV]');
    
    V(1) = Vfp(2)+0.1;

    for j=1:length(t)-1  
        k1v = (Iapp-Gl*(V(j)-El)-Gp*pinf(V(j))*(V(j)-Ena))/C;                 
        av = V(j)+k1v*dt;
        av = av+sqrt(2*D*dt)*eta(j);    
        k2v = (Iapp-Gl*(av-El)-Gp*pinf(av)*(av-Ena))/C;     
        V(j+1) = V(j)+(k1v+k2v)*dt/2;
        V(j+1) = V(j+1)+sqrt(2*D*dt)*eta(j);
    end
    
    V3 = V;
    
    figure(1)
    subplot(1,2,2)    
    plot(t,V,'-b','linewidth',2);
    plot(t,Vfp(2)*ones(1,length(t)),'--b','linewidth',2);
    
end

SMPFIG = 1;
if SMPFIG == 1
 

    % Sample figures

    close(1)

    figure(1001)
    hold on
    plot(vv,F(vv),'r','linewidth',2);
    plot([-100 100],[0 0],'--','Color',[.7 .7 .7],'linewidth',1);
    axis([-80 80 -10 30]);
    set(gca,'fontsize',20);
    xlabel('V  [mV]');
    ylabel('dV / dt');
    for k=1:length(Vfp)
        if Eig(k)<0
             plot(Vfp(k),0,'ok','linewidth',2);
             plot(Vfp(k),0,'xk','linewidth',2);
        elseif Eig(k)>0       
            plot(Vfp(k),0,'ok','linewidth',2);        
         end
    end

    if length(Vfp)==1     
        figure(1002)
        hold on
        plot(t,V,'-b','linewidth',2);
        axis([0 Tmax -80 40]);
        set(gca,'fontsize',24);
        xlabel('t  [ms]');
        ylabel('V  [mV]');
    elseif length(Vfp)==3
        figure(1002)
        hold on
        plot(t,V1,'-b','linewidth',2);
        plot(t,V3,'-b','linewidth',2);
        plot(t,Vfp(2)*ones(1,length(t)),'--b','linewidth',2);
        axis([0 Tmax -80 40]);
        set(gca,'fontsize',24);
        xlabel('t  [ms]');
        ylabel('V  [mV]');

    end
end

NSE = 0;
if NSE == 1
    
    
    D  = 5;
    
    Tmax = 100000;
    dt = 0.1;
    t = 0:dt:Tmax;
    eta = randn(1,length(t));

    
    V(1) = Vfp(2);

    for j=1:length(t)-1  
        k1v = (Iapp-Gl*(V(j)-El)-Gp*pinf(V(j))*(V(j)-Ena))/C;                 
        av = V(j)+k1v*dt;
        av = av+sqrt(2*D*dt)*eta(j);    
        k2v = (Iapp-Gl*(av-El)-Gp*pinf(av)*(av-Ena))/C;     
        V(j+1) = V(j)+(k1v+k2v)*dt/2;
        V(j+1) = V(j+1)+sqrt(2*D*dt)*eta(j);
    end
    
    figure
    hold on
    plot(t,V,'-b','linewidth',2);
    axis([0 Tmax -80 40]);
    set(gca,'fontsize',24);
    xlabel('t  [ms]');
    ylabel('V  [mV]');
end
