%
%  cncafib4ex.m   Steve Cox, 2/1/07
%
% solve the active fiber via sparse hybrid trapezoid
%
% with 3 calcium channels AND NaCa exchange
%
%  and solution of the u=[c;b] reaction diffusion system
%
%  and gKCa
%
%  and ER with Ryanodine Receptors
%
% usage:	[ca, t] = cncafib4ex(a,ell,N,dt,t1,t2,T,Iapp,pinc)
%
%		a = fiber radius (cm)
%		ell = fiber length (cm)
%       N = number of compartments
%       dt = timestep (ms)
%		t1 = start of current pulse (ms)
%		t2 = end of current pulse (ms)
%       Tfin = stopping time (ms)
%		Iapp = size of current pulse (micro amps)
%       pinc = number of time steps between plots
%
% e.g.:	[v, t] = cncafib3ex(.0001,.1,200,.02,2,4,15,.0002,10,0)
%

function [vmid, t] = cncafib3ex(rad,ell,N,dt,t1,t2,Tfin,Iapp,pinc)

Cm = 1;		% micro F / cm^2
R2 = 0.2; % 0.034;		% k Ohm cm
dx = ell/N;		% patch length
A = 2*pi*rad*dx;		% patch surface area
x = dx/2:dx:ell-dx/2;	% vector of patch midpoints
Nmid = round(N/2);

Nt = ceil(Tfin/dt)+1;
ICaT = zeros(Nt+1,1); ICaN = ICaT; ICaL = ICaT; vmid = ICaT;
ICa = ICaT; ca = ICa; bca = ICa; INaCa = ca;

E = struct('K', -77, 'Na', 56, 'Cl', -68);
g = struct('K', 36, 'Na', 120, 'Cl', 1/60, ...
           'CaT', .25, 'CaN', 2.5, 'CaL', 2.5, 'KCa', 10);
        
e1 = zeros(N,1);
e1(Nmid/2) = 1;
Iapp = Iapp*e1/A;
 
e = ones(N,1);

Na3 = (9/152)^3;
knaca = 4e2;   % micro Amp / cm^2

cr = .05*e;   % micro molar
co = 1e3;  % micro molar
BT = 5e2;    % micro molar
k1 = 1.5e-3;     % 1 / micro molar / ms
k2 = 0.3e-3;     % 1 / ms
br = BT*cr./(cr + k2/k1);
Dc = 220e-11;     % cm^2 / ms   parker
Db = 110e-11;     % cm^2 / ms    parker
F = 96485.3399;  % coulombs / mole   Faraday's constant
Jpmca = 1e-4;  % micro mole/s/cm^2
Jleak = Jpmca/21;
Kd = 1;    % micromolar

lam2 = rad/2/R2/g.Cl;

S = spdiags([-e 2*e -e], -1:1, N, N)/dx/dx;
S(1,1) = 1/dx/dx;
S(N,N) = 1/dx/dx;

eN = speye(N);
R = [-Dc*S-k1*BT*eN  k2*eN; k1*BT*eN -Db*S-k2*eN];
Rp = 2*speye(2*N) + dt*R;
Rm = 2*speye(2*N) - dt*R;
[LRm,URm] = lu(Rm);

B = lam2*g.Cl*S;

ci = cr;
Vr = fsolve(@(V) Iss(V,E,g,B,ci,co,knaca,Na3),-70*e);        % initial conditions

n1 = an(Vr)./(an(Vr)+bn(Vr)); n14 = n1.^4;
m1 = am(Vr)./(am(Vr)+bm(Vr)); 
h1 = ah(Vr)./(ah(Vr)+bh(Vr)); m13h1 = m1.^3.*h1;

m1T = amT(Vr)./(amT(Vr)+bmT(Vr)); 
h1T = ahT(Vr)./(ahT(Vr)+bhT(Vr)); m12h1T = m1T.^2.*h1T;

m1N = amN(Vr)./(amN(Vr)+bmN(Vr)); 
h1N = ahN(Vr)./(ahN(Vr)+bhN(Vr)); m12h1N = m1N.^2.*h1N;

m1L = amL(Vr)./(amL(Vr)+bmL(Vr)); m12L = m1L.^2;

inaca = -knaca*(Na3*exp(0.013*Vr)-(ci/co).*exp(-2*0.013*Vr));

m1KCa = amKCa(Vr,cr)./(amKCa(Vr,cr)+bmKCa(Vr,cr));

ICaT(1) = g.CaT*m12h1T(Nmid)*Phi(cr(Nmid),co,Vr(Nmid));
ICaN(1) = g.CaN*m12h1N(Nmid)*Phi(cr(Nmid),co,Vr(Nmid));
ICaL(1) = g.CaL*m12L(Nmid)*Phi(cr(Nmid),co,Vr(Nmid));
INaCa(1) = inaca(Nmid);
ICa(1) = ICaT(1)+ICaN(1)+ICaL(1)+INaCa(1);
IKCa(1) = g.KCa*m1KCa(Nmid)*(Vr(Nmid)-E.K);
ca(1) = cr(Nmid);
bca(1) = br(Nmid);
vmid(1) = Vr(Nmid);

u = [cr; br];           % initially at rest

ICa1 = (g.CaT*m12h1T + g.CaN*m12h1N + g.CaL*m12L).*Phi(u(1:N),co,Vr);

v = Vr;
a = an(v);  b = bn(v);
n2 = ( (2/dt-a-b).*n1 + 2*a) ./ (2/dt + a + b); n24 = n2.^4;

a = am(v);  b = bm(v);
m2 = ( (2/dt-a-b).*m1 + 2*a) ./ (2/dt + a + b);

a = ah(v);  b = bh(v);
h2 = ( (2/dt-a-b).*h1 + 2*a) ./ (2/dt + a + b); m23h2 = m2.^3.*h2;

a = amT(v);  b = bmT(v);
m2T = ( (2/dt-a-b).*m1T + 2*a) ./ (2/dt + a + b);
a = ahT(v);  b = bhT(v);
h2T = ( (2/dt-a-b).*h1T + 2*a) ./ (2/dt + a + b); m22h2T = m2T.^2.*h2T;

a = amN(v);  b = bmN(v);
m2N = ( (2/dt-a-b).*m1N + 2*a) ./ (2/dt + a + b);
a = ahN(v);  b = bhN(v);
h2N = ( (2/dt-a-b).*h1N + 2*a) ./ (2/dt + a + b); m22h2N = m2N.^2.*h2N;

a = amL(v);  b = bmL(v);
m2L = ( (2/dt-a-b).*m1L + 2*a) ./ (2/dt + a + b); m22L = m2L.^2;

a = amKCa(v,cr);  b = bmKCa(v,cr);
m2KCa = ( (2/dt-a-b).*m1KCa + 2*a) ./ (2/dt + a + b);

inaca = -knaca*(Na3*exp(0.013*v)-(u(1:N)/co).*exp(-2*0.013*v));

ICaT(2) = g.CaT*m22h2T(Nmid)*Phi(u(Nmid),co,Vr(Nmid));
ICaN(2) = g.CaN*m22h2N(Nmid)*Phi(u(Nmid),co,Vr(Nmid));
ICaL(2) = g.CaL*m22L(Nmid)*Phi(u(Nmid),co,Vr(Nmid));
IKCa(2) = g.KCa*m2KCa(Nmid)*(Vr(Nmid)-E.K);
INaCa(2) = inaca(Nmid);
ICa(2) = ICaT(2)+ICaN(2)+ICaL(2)+INaCa(2);
ca(2) = u(Nmid);
bca(2) = u(Nmid+N);
vmid(2) = Vr(Nmid);

t = 0;
if pinc > 0
   figure(1)
   plot3(x,t*ones(N,1),u(N+1:2*N))
   hold on
   figure(2)
   plot3(x,t*ones(N,1),u(1:N))
   hold on
   figure(3)
   plot3(x,t*ones(N,1),-ICa1)
   hold on
   figure(4)
   plot3(x,t*ones(N,1),v)
   hold on
end

i1 = 0;
t = dt;
i2 = Iapp*(t>t1)*(t<t2);

dB = diag(B);

tf = 2*Cm/dt;

d2 = g.Na*m23h2 + g.K*n24 + g.Cl + g.KCa*m2KCa;

ICa2 = (g.CaT*m22h2T + g.CaN*m22h2N + g.CaL*m22L).*Phi(u(1:N),co,v);
   
f = g.Na*(m23h2 + m13h1)*E.Na + g.K*(n24+n14)*E.K + 2*g.Cl*E.Cl ...
      - 2*inaca + i2 + i1 - (ICa1 + ICa2) + g.KCa*(m2KCa+m1KCa)*E.K;

r = (tf - d2).*v - B*v + f;

B(1:N+1:end) = dB + d2 + tf;         % update the diagonal

v = B\r;   % v at dt

if pinc > 0
   figure(1)
   plot3(x,t*ones(N,1),u(N+1:2*N))
   figure(2)
   plot3(x,t*ones(N,1),u(1:N))
   figure(3)
   plot3(x,t*ones(N,1),-ICa2)
   figure(4)
   plot3(x,t*ones(N,1),v)
end

i1 = i2; n1 = n2; m1 = m2; h1 = h2; d1 = d2; m13h1 = m23h2; n14 = n24; 
h1N = h2N; m1N = m2N; h1T = h2T; m1T = m2T; m1L = m2L; ICa1 = ICa2;
m1KCa = m2KCa;

for j=2:Nt,
    
    p1 = -(2/rad)*(ICa2+inaca)/(2*F);
    p2 = (2/rad)*(Jpmca*u(1:N)./(Kd+u(1:N))-Jleak);
    p = p1 - p2;
    tmp = k1*u(1:N).*u(N+1:2*N);
    q = 2*[tmp + p; -tmp];
    u = URm \ ( LRm \ (Rp*u + dt*q) );
    
    t = j*dt;

    i2 = Iapp*(t>t1)*(t<t2);

    a = an(v);  b = bn(v);
    n2 = ( (2/dt-a-b).*n1 + 2*a) ./ (2/dt + a + b); n24 = n2.^4;

    a = am(v);  b = bm(v);
    m2 = ( (2/dt-a-b).*m1 + 2*a) ./ (2/dt + a + b);    

    a = ah(v);  b = bh(v);
    h2 = ( (2/dt-a-b).*h1 + 2*a) ./ (2/dt + a + b); m23h2 = m2.^3.*h2;

    a = amT(v);  b = bmT(v);
    m2T = ( (2/dt-a-b).*m1T + 2*a) ./ (2/dt + a + b);
    a = ahT(v);  b = bhT(v);
    h2T = ( (2/dt-a-b).*h1T + 2*a) ./ (2/dt + a + b); m22h2T = m2T.^2.*h2T;

    a = amN(v);  b = bmN(v);
    m2N = ( (2/dt-a-b).*m1N + 2*a) ./ (2/dt + a + b);
    a = ahN(v);  b = bhN(v);
    h2N = ( (2/dt-a-b).*h1N + 2*a) ./ (2/dt + a + b); m22h2N = m2N.^2.*h2N;

    a = amL(v);  b = bmL(v);
    m2L = ( (2/dt-a-b).*m1L + 2*a) ./ (2/dt + a + b); m22L = m2L.^2;
    
    a = amKCa(v,u(1:N));  b = bmKCa(v,u(1:N));
    m2KCa = ( (2/dt-a-b).*m1KCa + 2*a) ./ (2/dt + a + b);

    d2 = g.Na*m23h2 + g.K*n24 + g.Cl + g.KCa*m2KCa;

    ICa2 = (g.CaT*m22h2T + g.CaN*m22h2N + g.CaL*m22L).*Phi(u(1:N),co,v);
    
    inaca = -knaca*(Na3*exp(0.013*v)-(u(1:N)/co).*exp(-2*0.013*v));

    f = g.Na*(m23h2 + m13h1)*E.Na + g.K*(n24+n14)*E.K + 2*g.Cl*E.Cl ...
        - 2*inaca + i2 + i1 - (ICa1 + ICa2) + g.KCa*(m2KCa+m1KCa)*E.K;

    r = (2*tf - d2 + d1).*v - r + f;

    B(1:N+1:end) = dB + d2 + tf;         % update the diagonal

    v = B\r;

    if mod(j,pinc) == 0
        figure(1)
        plot3(x,t*ones(N,1),u(N+1:2*N))
        figure(2)
        plot3(x,t*ones(N,1),u(1:N))
        figure(3)
        plot3(x,t*ones(N,1),-ICa2)
        figure(4)
        plot3(x,t*ones(N,1),v)
        disp([max(abs(p1)) max(abs(p2))])
    end

    ICaT(j+1) = g.CaT*m22h2T(Nmid)*Phi(u(Nmid),co,v(Nmid));
    ICaN(j+1) = g.CaN*m22h2N(Nmid)*Phi(u(Nmid),co,v(Nmid));
    ICaL(j+1) = g.CaL*m22L(Nmid)*Phi(u(Nmid),co,v(Nmid));
    IKCa(j+1) = g.KCa*m2KCa(Nmid)*(v(Nmid)-E.K);
    INaCa(j+1) = inaca(Nmid);
    vmid(j+1) = v(Nmid);
    ICa(j+1) = ICaT(j+1)+ICaN(j+1)+ICaL(j+1)+INaCa(j+1);
    ca(j+1) = u(Nmid);
    bca(j+1) = u(Nmid+N);

    i1 = i2; n1 = n2; m1 = m2; h1 = h2; d1 = d2; m13h1 = m23h2; n14 = n24;
    h1N = h2N; m1N = m2N; h1T = h2T; m1T = m2T; m1L = m2L; ICa1 = ICa2;
    m1KCa = m2KCa;
    
end

if pinc > 0
    figure(1)
   xlabel('x  (cm)','fontsize',16)
   ylabel('t  (ms)','fontsize',16)
   zlabel('b  (\muM)','fontsize',16)
   hold off
   figure(2)
   xlabel('x  (cm)','fontsize',16)
   ylabel('t  (ms)','fontsize',16)
   zlabel('c  (\muM)','fontsize',16)
   hold off
   figure(3)
   xlabel('x  (cm)','fontsize',16)
   ylabel('t  (ms)','fontsize',16)
   zlabel('-I_{Ca}  (\muA/cm^2)','fontsize',16)
   hold off
   figure(4)
   xlabel('x  (cm)','fontsize',16)
   ylabel('t  (ms)','fontsize',16)
   zlabel('v  (mV)','fontsize',16)
   hold off
end

t = linspace(0,Tfin,length(ICaT))';

figure(2)
subplot(4,1,1)
plot(t,vmid,'linewidth',2)
ylabel('v (mV)','fontsize',12)
%ylim([-80 60])

subplot(4,1,2)
plot(t,IKCa,'linewidth',2)
ylabel('I_{K,Ca} (\muA/cm^2)','fontsize',12)
%ylim([-75 10])

subplot(4,1,3)
plot(t,ca,'linewidth',2)
ylabel('c (\muM)','fontsize',12)
%ylim([0 45])

subplot(4,1,4)
plot(t,bca,'linewidth',2)
xlabel('t  (ms)','fontsize',12)
ylabel('b (\muM)','fontsize',12)

return

function val = Iss(V,E,g,B,ci,co,knaca,Na3)
val = B*V + g.Na*(am(V)./(am(V)+bm(V))).^3.*(ah(V)./(ah(V)+bh(V))).*(V-E.Na) + ...
            g.K*(an(V)./(an(V)+bn(V))).^4.*(V-E.K) + g.Cl.*(V-E.Cl) + ...
   (g.CaT*(amT(V)./(amT(V)+bmT(V))).^2.*(ahT(V)./(ahT(V)+bhT(V))) + ...
    g.CaN*(amN(V)./(amN(V)+bmN(V))).^2.*(ahN(V)./(ahN(V)+bhN(V))) + ...
    g.CaL*(amL(V)./(amL(V)+bmL(V))).^2) .*Phi(ci,co,V)- ...
    knaca*(Na3*exp(0.013*V)-(ci/co).*exp(-2*0.013*V)) + ...
    g.KCa*(amKCa(V,ci)./(amKCa(V,ci)+bmKCa(V,ci))).*(V-E.K);

function val = an(v)
val = .01*(10-(v+71))./(exp(1-(v+71)/10)-1);

function val = bn(v)
val = .125*exp(-(v+71)/80);

function val = am(v)
val = .1*(25-(v+71))./(exp(2.5-(v+71)/10)-1);

function val = bm(v)
val = 4*exp(-(v+71)/18);

function val = ah(v)
val = 0.07*exp(-(v+71)/20);

function val = bh(v)
val = 1./(exp(3-(v+71)/10)+1);

function val = Phi(ci,co,v)
F = 96485; % C/mol
R = 8.314;  % J/deg K/mole
T = 300;   % deg K
e = exp(2*F*v*1e-3/R/T);
val = v.*(1-(ci/co).*e)./(1-e);

function val = amL(V)
a_L = 15.69; b_L = 81.5; c_L = 0.29; d_L = 10.86;
val = a_L*(b_L-V)./(exp((b_L-V)/10)-1);

function val = bmL(V)
a_L = 15.69; b_L = 81.5; c_L = 0.29; d_L = 10.86;
val = c_L*exp(-V/d_L);

function val = amN(V)
a_N = 0.19; b_N = 19.88; c_N = 0.046; d_N = 20.73;
e_N = 1.6e-4; f_N=48.46; g_N=39;
val = a_N*(b_N-V)./(exp((b_N-V)/10)-1);

function val = bmN(V)
a_N = 0.19; b_N = 19.88; c_N = 0.046; d_N = 20.73;
e_N = 1.6e-4; f_N=48.46; g_N=39;
val = c_N*exp(-V/d_N);

function val = ahN(V)
a_N = 0.19; b_N = 19.88; c_N = 0.046; d_N = 20.73;
e_N = 1.6e-4; f_N=48.46; g_N=39;
val = e_N*exp(-V/f_N);

function val = bhN(V)
a_N = 0.19; b_N = 19.88; c_N = 0.046; d_N = 20.73;
e_N = 1.6e-4; f_N=48.46; g_N=39;
val  = 1./(1+exp((g_N-V)/10));

function val = amT(V)
a_T = 0.2; b_T = 19.26; c_T = 0.009; d_T = 22.03;
e_T = 1e-6; f_T=16.26; g_T=29.79;
val = a_T*(b_T-V)./(exp((b_T-V)/10)-1);

function val = bmT(V)
a_T = 0.2; b_T = 19.26; c_T = 0.009; d_T = 22.03;
e_T = 1e-6; f_T=16.26; g_T=29.79;
val = c_T*exp(-V/d_T);

function val = ahT(V)
a_T = 0.2; b_T = 19.26; c_T = 0.009; d_T = 22.03;
e_T = 1e-6; f_T=16.26; g_T=29.79;
val = e_T*exp(-V/f_T);

function val = bhT(V)
a_T = 0.2; b_T = 19.26; c_T = 0.009; d_T = 22.03;
e_T = 1e-6; f_T=16.26; g_T=29.79;
val = 1./(1+exp((g_T-V)/10));

function val = amKCa(v,c)
FoRT = 96.485/8.314/300; 
val = 0.28*c./(c+.48*exp(-2*.84*v*FoRT));

function val = bmKCa(v,c)
FoRT = 96.485/8.314/300;
val = 0.48./(1+(c/.13e-3).*exp(2*v*FoRT));