% Author: David Catherall; Duke University
% Created: September 2017
% Description: Matlab Voltage Clamp Implementation for Schild 1994 for
% primary current, using constant calcium ion concentrations

clear
close all
%% Constants

% Schild 1994 Voltage Clamps Figure Values
% T=296;
% shiftca=0;
% shiftnaf=0;
% shiftnas=0;
% shiftk=0;
% gCat=2.75e-9;
% gNaf=29.5e-9;
% gNas=25.5e-9;
% gKd=30.0e-9;
% gKa=28.0e-9;
% gKds=11.5e-9;
% gCan=22.5e-9;

% C-Fiber Values
T=311;
shiftca=-7;
shiftnaf=-17.5;
shiftnas=-20.0;
shiftk=3.0;
gNaf=1.95e-6;
gNas=0.0295e-6;
gCat=0.00035e-6;
gCan=0.003e-6;
gKd=0.0051e-6;
gKa=0.004e-6;
gKds=0.003e-6;

R=8314;
F=96500;
Nai=8.9;
Nao=154;
Ena=R*T/(F)*log(Nao/Nai);
Ko=5.4;
Ki=145.0;
Ek=R*T/(F)*log(Ko/Ki);
Cai=0.000117;
Cao=2.0;
Eca=R*T/(2*F)*log(Cao/Cai)-78.7;

%% Cat
V=-60:10:-20;
%Analytical
for j=1:length(V)
    Q10dt=1.90;
    Q10ft=2.20;
    dt0=1.0/(1.0+exp((-80+54.00+shiftca)/-5.75));
    ft0=1.0/(1.0+exp((-80+68.00+shiftca)/6.0));
    dtinf=1.0/(1.0+exp((V(j)+54.00+shiftca)/-5.75));
    taudt=22*exp(-(0.052)^2*(V(j)+68.0)^2)+2.5;
    taudt=taudt*Q10dt^((22-T+274)/10);
    ftinf=1.0/(1.0+exp((V(j)+68.00+shiftca)/6.0));
    tauft=103*exp(-(0.050)^2*(V(j)+58.0)^2)+12.5;
    tauft=tauft*Q10ft^((22-T+274)/10);

    t=0:0.005:200-0.005;

    dt=(dt0-dtinf)*exp(-t/taudt)+dtinf;
    ft=(ft0-ftinf)*exp(-t/tauft)+ftinf;

    Icat(j,:)=gCat.*dt.*ft.*(V(j)-Eca)*10^9;

    dtime=t(2);
    dtnum(1)=dt0;
    ftnum(1)=ft0;
    ddtnum(1)=0;
    dftnum(1)=0;
% Numerical
    for i=2:length(t)
        ddtnum(i)=(dtinf-dtnum(i-1))/taudt;
        dftnum(i)=(ftinf-ftnum(i-1))/tauft;
        dtnum(i)=dtnum(i-1)+ddtnum(i)*dtime;
        ftnum(i)=ftnum(i-1)+dftnum(i)*dtime;
    end

    Icatnum(j,:)=gCat.*dtnum.*ftnum.*(V(j)-Eca)*10^9;
end

%% Naf

V=[-40:10:0];
% Analytical
for j=1:length(V)
    Q10m=2.30;
    Q10h=1.50;
    mf0=1.0/(1.0+exp((-80+41.35+shiftnaf)/-4.75));
    hf0=1.0/(1.0+exp((-80+62.00+shiftnaf)/4.5));
    j0=1.0/(1.0+exp((-80+40.00)/1.50));
    mfinf=1.0/(1.0+exp((V(j)+41.35+shiftnaf)/-4.75));
    taumf=0.75*exp(-(0.0635)^2*(V(j)+40.35)^2)+0.12;
    taumf=taumf*Q10m^((22-T+274)/10);
    hfinf=1.0/(1.0+exp((V(j)+62.00+shiftnaf)/4.5));
    tauhf=6.5*exp(-(0.0295)^2*(V(j)+75.00)^2)+0.55;
    tauhf=tauhf*Q10h^((22-T+274)/10);
    jinf=1.0/(1.0+exp((V(j)+40.00)/1.50));
    tauj=(25.0/(1.0+exp((V(j)-20.00)/4.50)))+0.01;

    t=0:0.005:10-0.005;

    mf=(mf0-mfinf)*exp(-t/taumf)+mfinf;
    hf=(hf0-hfinf)*exp(-t/tauhf)+hfinf;
    jf=(j0-jinf)*exp(-t/tauj)+jinf;

    Inaf(j,:)=gNaf.*mf.^3.*hf.*jf.*(V(j)-Ena)*10^9;

    dtime=t(2);
    mfnum(1)=mf0;
    hfnum(1)=hf0;
    jnum(1)=j0;
    dmfnum(1)=0;
    dhfnum(1)=0;
    djnum(1)=0;
% Numerical
    for i=2:length(t)
        dmfnum(i)=(mfinf-mfnum(i-1))/taumf;
        dhfnum(i)=(hfinf-hfnum(i-1))/tauhf;
        djnum(i)=(jinf-jnum(i-1))/tauj;
        mfnum(i)=mfnum(i-1)+dmfnum(i)*dtime;
        hfnum(i)=hfnum(i-1)+dhfnum(i)*dtime;
        jnum(i)=jnum(i-1)+djnum(i)*dtime;
    end

    Inafnum(j,:)=gNaf.*mfnum.^3.*hfnum.*jnum.*(V(j)-Ena)*10^9;
end

% Nas

V=-10:10:30;
% Analytical
for j=1:length(V)
    ms0=1.0/(1.0+exp((-80+20.35+shiftnas)/-4.45));
    hs0=1.0/(1.0+exp((-80+18.00+shiftnas)/4.5));
    msinf=1.0/(1.0+exp((V(j)+20.35+shiftnas)/-4.45));
    taums=1.50*exp(-(0.0595)^2*(V(j)+20.35)^2)+0.15;
    taums=taums*Q10m^((22-T+274)/10);
    hsinf=1.0/(1.0+exp((V(j)+18.00+shiftnas)/4.5));
    tauhs=4.95*exp(-(0.0335)^2*(V(j)+20.00)^2)+0.75;
    tauhs=tauhs*Q10h^((22-T+274)/10);

    t=0:0.005:10-0.005;

    ms=(ms0-msinf)*exp(-t/taums)+msinf;
    hs=(hs0-hsinf)*exp(-t/tauhs)+hsinf;

    Inas(j,:)=gNas.*ms.^3.*hs.*(V(j)-Ena)*10^9;

    dtime=t(2);
    msnum(1)=ms0;
    hsnum(1)=hs0;
    dmsnum(1)=0;
    dhsnum(1)=0;
% Numerical
    for i=2:length(t)
        dmsnum(i)=(msinf-msnum(i-1))/taums;
        dhsnum(i)=(hsinf-hsnum(i-1))/tauhs;
        msnum(i)=msnum(i-1)+dmsnum(i)*dtime;
        hsnum(i)=hsnum(i-1)+dhsnum(i)*dtime;
    end

    Inasnum(j,:)=gNas.*msnum.^3.*hsnum.*(V(j)-Ena)*10^9;
end

% Kd

V=-40:10:40;
% Analytical
for j=1:length(V)
    Q10n=1.40;
    n0=1.0/(1.0+exp((-80+14.62+shiftk)/-18.38));
    ninf=1.0/(1.0+exp((V(j)+14.62+shiftk)/-18.38));
    an=(0.001265*(V(j)+14.273))/(1.0-exp((V(j)+14.273)/-10.0));
    bn=0.125*exp((V(j)+55.0)/-2.5);
    taun=(1/(an+bn))+1;
    taun=taun*Q10n^((22-T+274)/10);

    t=0:0.005:350-0.005;

    n=(n0-ninf)*exp(-t/taun)+ninf;

    Ikd(j,:)=gKd.*n.*(V(j)-Ek)*10^9;

    dtime=t(2);
    nnum(1)=n0;
    dnnum(1)=0;
% Numerical
    for i=2:length(t)
        dnnum(i)=(ninf-nnum(i-1))/taun;
        nnum(i)=nnum(i-1)+dnnum(i)*dtime;
    end

    Ikdnum(j,:)=gKd.*nnum.*(V(j)-Ek)*10^9;
end

% Ktrans

V=-40:10:30;
% Analytical 
for j=1:length(V)
    Q10kt=1.93;
    p0=1.0/(1.0+exp((-80+28.0+shiftk)/-28.0));
    q0=1.0/(1.0+exp((-80+58.0+shiftk)/7.0));
    pinf=1.0/(1.0+exp((V(j)+28.0+shiftk)/-28.0));
    taup=5.0*exp(-(0.022)^2*(V(j)+65.0)^2)+2.5;
    taup=taup*Q10kt^((22-T+274)/10);
    qinf=1.0/(1.0+exp((V(j)+58.0+shiftk)/7.0));
    tauq=100.0*exp(-(0.035)^2*(V(j)+30.00)^2)+10.5;
    tauq=tauq*Q10kt^((22-T+274)/10);

    t=0:0.005:600-0.005;

    p=(p0-pinf)*exp(-t/taup)+pinf;
    q=(q0-qinf)*exp(-t/tauq)+qinf;

    Ika(j,:)=gKa.*p.^3.*q.*(V(j)-Ek)*10^9;

    dtime=t(2);
    pnum(1)=p0;
    qnum(1)=q0;
    dpnum(1)=0;
    dqnum(1)=0;
% Numerical
    for i=2:length(t)
        dpnum(i)=(pinf-pnum(i-1))/taup;
        dqnum(i)=(qinf-qnum(i-1))/tauq;
        pnum(i)=pnum(i-1)+dpnum(i)*dtime;
        qnum(i)=qnum(i-1)+dqnum(i)*dtime;
    end

    Ikanum(j,:)=gKa.*pnum.^3.*qnum.*(V(j)-Ek)*10^9;
end
% Analytical
for j=1:length(V)
    x0=1.0/(1.0+exp((-80+39.59+shiftk)/-14.68));
    y0=1.0/(1.0+exp((-80+48.0+shiftk)/7.0));
    xinf=1.0/(1.0+exp((V(j)+39.59+shiftk)/-14.68));
    taux=5.0*exp(-(0.022)^2*(V(j)+65.0)^2)+2.5;
    taux=taux*Q10kt^((22-T+274)/10);
    yinf=1.0/(1.0+exp((V(j)+48.0+shiftk)/7.0));
    tauy=7500.0;
    tauy=tauy*Q10kt^((22-T+274)/10);

    t=0:0.005:600-0.005;

    x=(x0-xinf)*exp(-t/taux)+xinf;
    y=(y0-yinf)*exp(-t/tauy)+yinf;

    Ikds(j,:)=gKds.*x.^3.*y.*(V(j)-Ek)*10^9;

    dtime=t(2);
    xnum(1)=x0;
    ynum(1)=y0;
    dxnum(1)=0;
    dynum(1)=0;
% Numerical
    for i=2:length(t)
        dxnum(i)=(xinf-xnum(i-1))/taux;
        dynum(i)=(yinf-ynum(i-1))/tauy;
        xnum(i)=xnum(i-1)+dxnum(i)*dtime;
        ynum(i)=ynum(i-1)+dynum(i)*dtime;
    end

    Ikdsnum(j,:)=gKds.*xnum.^3.*ynum.*(V(j)-Ek)*10^9;
end

Iktrans=Ika+Ikds;
Iktransnum=Ikanum+Ikdsnum;

% Can

V=-20:10:20;
% Analytical
for j=1:length(V)
    Q10can=4.30;
    dn0=1.0/(1.0+exp((-80+20.0+shiftca)/-4.5));
    fn10=1.0/(1.0+exp((-80+20.0+shiftca)/25.0));
    rn0=0.2/(1.0+exp((-80+5.0+shiftca)/-10.0));
    fn20=rn0+1.0/(1.0+exp((-80+40.0+shiftca)/10.0));
    dninf=1.0/(1.0+exp((V(j)+20.0+shiftca)/-4.5));
    taudn=3.25*exp(-(0.042)^2*(V(j)+31.0)^2)+0.395;
    taudn=taudn*Q10can^((22-T+274)/10);
    fn1inf=1.0/(1.0+exp((V(j)+20.0+shiftca)/25.0));
    taufn1=33.5*exp(-(0.0395)^2*(V(j)+30.0)^2)+5.0;
    taufn1=taufn1*Q10can^((22-T+274)/10);
    rn=0.2/(1.0+exp((V(j)+5.0+shiftca)/-10.0));
    fn2inf=rn+1.0/(1.0+exp((V(j)+40.0+shiftca)/10.0));
    taufn2=225.0*exp(-(0.0275)^2*(V(j)+40.0)^2)+75.0;
    taufn2=taufn2*Q10can^((22-T+274)/10);

    t=0:0.005:200-0.005;

    dn=(dn0-dninf)*exp(-t/taudn)+dninf;
    fn1=(fn10-fn1inf)*exp(-t/taufn1)+fn1inf;
    fn2=(fn20-fn2inf)*exp(-t/taufn2)+fn2inf;

    Ican(j,:)=gCan.*dn.*(0.55*fn1+0.45*fn2).*(V(j)-Eca)*10^9;

    dtime=t(2);
    ddnnum(1)=dn0;
    fn1num(1)=fn10;
    fn2num(1)=fn20;
    ddnnum(1)=0;
    dfn1num(1)=0;
    dfn2num(1)=0;
% Numerical
    for i=2:length(t)
        ddnnum(i)=(dninf-ddnnum(i-1))/taudn;
        dfn1num(i)=(fn1inf-fn1num(i-1))/taufn1;
        dfn2num(i)=(fn2inf-fn2num(i-1))/taufn2;
        ddnnum(i)=ddnnum(i-1)+ddnnum(i)*dtime;
        fn1num(i)=fn1num(i-1)+dfn1num(i)*dtime;
        fn2num(i)=fn2num(i-1)+dfn2num(i)*dtime;
    end

    Icannum(j,:)=gCan.*ddnnum.*(0.55*fn1num+0.45*fn2num).*(V(j)-Eca)*10^9;
end

% save('33CMatlabClamp/ConstClamps','Ican','Icat','Ikd','Iktrans','Inaf','Inas')