% TBI (traumatic brain injury block for bistable equation
%This code calculates the critical relationship for A_r vs Y (the
%  length for propagation block following TBI)
% 
% here we use the HH upstroke function
function block
global  Ar un F_int U1 Rm gnabar
 
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

 
dv = 0.01;
V0 = [-12:dv:115];
gnabar = 120;  %nominal value
gnabar = 120;
% first, estimate the slope of I_ion at the rest state

%find 3 zeros
Ion_current = HH_rhs(V0);
%ndx has three numbers
ndx = find(Ion_current(1:end-1).*Ion_current(2:end)<=0);
   
 % slope at the rest point (for proper scaling)
 RM_inv = (Ion_current(ndx(1)+1)-Ion_current(ndx(1)))/dv;
 % this is the membrane resisitance
 Rm = 1/RM_inv
 
%%Specify (or determine) the three zeros of f(u)

%use bisection to determine the three zeros
a = V0(ndx);
b = V0(ndx+1);
r = bisect(a,b,@HH_rhs);

U0 = r(1);
U1=r(2);
U2 = r(3);

N = 1000; % discretization of the u interval
du = (U2-U0)/N;
un = [U0:du:U2];

% calculate the integral of f
 F_int = cumsum(feval(un))*du;

 %find the zero of F_int using bisection
 % guess a lower bound
 a = U0+du;
 % guess an upper bound
 b=U2;
 
umx = bisect(a,b,@fint);  %This is the zero of F_int
  
%plot the two basic trajectories
Fa = fint(U1);
u1=un(2:end);
F1 = fint(U2);
w1 = sqrt(2*(F1-fint(u1)));
 w1s=sqrt(2*(F1-Fa));
 % discretize the (U0,umx) interval
  
%The two basic trajectories
Fa = fint(U1);
% the upper trajectory
u1=un;
F1 = fint(U2);
w1 = sqrt(2*(F1-fint(u1)));
 w1s=sqrt(2*(F1-Fa));
 
  % discretize the (U0,umx) interval
  %the lower trqjectory
  M = 50; %
 du0 = (umx-U0)/M;
u0 = [U0:du0:umx];
 Fw = fint(u0);
w0 = sqrt(-2*Fw);
 
MN = 100;
% discretize the (U0,U1) interval for starting values Us
dua = (U1-U0)/MN;
F1 = fint(U2);
Fa = fint(U1);
Uj = [U0+dua:dua:U1-6*dua];  %These are the starting values for the intermediate trajectory
for j = 1:length(Uj)
    Us = Uj(j);
    
 % find the starting value for W
Ws0 = sqrt(-2*fint(Us));

Ar =(F1-Fa-Ws0^2/2)/(fint(Us)-Fa); % the critical value of As
A(j) = Ar;

tspan=[0:1:20000];
s = [Us Ws0];
% integrate the intermediate curve
[T,S] = ode23s(@deRHS,tspan,s,odeset('Events',@events)); 

Y(j) = T(end);
end

[miny,minj]= min(Y)
minA = A(minj)
 
figure(3)
plot(Uj,Y,Uj(minj),Y(minj),'*')
xlabel('Us' )
xlabel( 'Y')

figure(4)
semilogy(Uj,A,Uj(minj),A(minj),'*')
xlabel('Us' )
xlabel( 'A_r')



figure(1)
semilogy(Y(1:minj),A(1:minj),Y(minj:end),A(minj:end),'--',miny,minA,'k*', ...
     'linewidth',2)
 
ylabel('A_r','fontsize',20)
xlabel('Y','fontsize',20)
  
% for the minimal Y value, make a phase portrait


% Pick a starting value for Us < U1
Us = Uj(minj);  
% find the corresponding starting value for W
Ws0 = sqrt(-2*fint(Us));
Ar =(F1-Fa-Ws0^2/2)/(fint(Us)-Fa);

Ks = Ws0^2/2+Ar*fint(Us);
 Fmsk = (Ks-Ar*fint(u0) >= 0);
  Ws = sqrt(2* (Ks-Ar*fint(u0(Fmsk))));
 
tspan=[0:1:10000];  % hopefully this time interval is sufficiently long
s = [Us Ws0];
% integrate the intermediate curve
[T,S] = ode23s(@deRHS,tspan,s,odeset('Events',@events)); 
figure(2)
semilogx(u0,w0,u1,w1 ,u0(Fmsk),Ws,'k--',umx,0,'k*', Us,Ws0,'k*',  ...
    U1,w1s,'k*',S(:,1),S(:,2),'g','linewidth',2)
xlabel('U','fontsize',20)
ylabel('W','fontsize',20)
%axis([-10  20 0 1400])
 

function s_prime=deRHS(t,s)  % this is the differential equation for the intermediate section

global  Ar
 u = s(1);
 w = s(2);
 Fu = w/Ar;
 Fw = - feval(u);
 s_prime = [Fu;Fw];
    
% here is where to specify the function f(u);
function f=feval(u)
global Rm
 
f = -Rm*HH_rhs(u);

        
function Ion_current = HH_rhs(V)
global gnabar
gkbar = 36;

h = 0.5961;
n = 0.3177;

AM=0.1*(25.-V)./(exp(0.1*(25.0-V))-1.);
BM=4.0*exp(-V/18);
m = AM./(AM+BM);

gna=gnabar*m.^3.*h;
 
gk=gkbar*n.^4;

Vna = 115 ;
Vk =  -12 ;
Vl =  10.5988;
Icl=0.3*(V-Vl);

Ion_current = -gna.*(V-Vna)-gk.*(V-Vk)-Icl;

function fI = fint(u)
global un F_int

%Calculate the integral of f
% may be a vector
for j = 1:length(u)
ndx = min(find(u(j)<=un));

if(ndx==1)
    fI(j) = 0;
else
%linear interpolation
u1 = un(ndx-1);
f1 = F_int(ndx-1);
u2 = un(ndx);
f2 = F_int(ndx);
fI(j) = ((u(j)-u1)*f2-(u(j)-u2)*f1)/(u2-u1);
end
end

% event function
function [value,isterminal,direction]=events(t,y)
global U1

value = y(1)-U1; 
isterminal = 1;
direction = 0;

% reminder how to use event functions
%[T,S] = ode23s(@deRHS,tspan,s,odeset('Events',@events)); 

function root = bisect(a,b,feval)

ul = a;
fl = feval(ul);
uu = b;
fu = feval(uu);

% we make the assumption, without checking, that 
% fu*fl<0

% if not, the algoritm fails to find a root.

N = 20;  % number of iterates
% the main bisection algorithm
for j = 1:N
uc = (ul+uu)/2;
fc = feval(uc);
ftest = (fc.*fl>=0);
ul = ftest.*uc+(1-ftest).*ul;
fl = ftest.*fc + (1-ftest).*fl;

uu = (1-ftest).*uc+ ftest.*uu;
fu = (1-ftest).*fc + ftest.*fu;
end
root = uc;