%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Trovato2020 model for human cardiac Purkinje electrophysiolgy           %
%                                                                         %
% Model description and validation in:                                    %
% C. Trovato, E. Passini, N. Nagy, et al., Human Purkinje in              %
% silico model enables mechanistic investigations into automaticity       %
% and pro-arrhythmic abnormalities, Journal of Molecular and Cellular     %
% Cardiology (2019), https://doi.org/10.1016/j.yjmcc.2020.04.001 .        %
% [open-access]                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dX=modTrovato2020(t,X) 
%% Extracellular Ionic Concentrations [mM]
cEx = [1.8, 140.0, 5.4];
cao = cEx(1);  %[Ca]o mM
nao = cEx(2);  %[Na]o mM
ko  = cEx(3);  %[K]o  mM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% State Variables
%Membrane Potential V
v = X(1);
%Ionic Intracellular Concentrations
nai = X(2);   nasl = X(3);   nass = X(4);  
ki = X(5);    ksl = X(6);    kss = X(7);      
cai = X(8);   casl = X(9);   cass = X(10);
cansr = X(11);   cajsr = X(12); cacsr=X(13);
%INa 
m = X(14);      hf  = X(15);   hs = X(16);
j = X(17);      hsp = X(18);   jp = X(19);
%INaL 
mL = X(20);     hL = X(21);    hLp = X(22); 
%Ito
a  = X(23); i=X(24);     i2=X(25);    
%ICaL
d    = X(26);   
ff   = X(27);  fs    = X(28); fcaf = X(29);   fcas = X(30);  
jca   = X(31); nca = X(32); ffp  = X(33);  fcafp = X(34); 
%ICaT
b=X(35);     g=X(36);
%IKr
xrf = X(37);    xrs = X(38);
%IKs
xs1 = X(39);    xs2 = X(40);
%If
y=X(41);
%IK1
xk1 = X(42);
%CaMKt
CaMKt=X(43);
%Releases      
u=X(44);    Jrel1=X(45);  Jrel2=X(46);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Physical Constants:
R = 8314.0;   % J/kmol/K
T = 310.0;    % K
F = 96485.0;  % C/mol
vffrt = v*F*F/(R*T);
vfrt  = v*F/(R*T);
%% Valence of ions
zna = 1;
zk  = 1;
zca = 2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Cell Geometry from ORd
% Cell geometry was approxymate by a cylinder of length L and radius r
L = 0.0164;                         % cm
rad = 0.00175;                      % cm
vcell = 1000*pi*rad^2*L;            % uL
% Geometric Area
Ageo = 2*pi*rad^2 + 2*pi*rad*L;     % cm^2
% Capacitive Area
Acap = 2*Ageo;                      % cm^2
% Compartment Volumes (6)
vmyo = 0.60*vcell;                  % uL
vnsr = 0.04*vcell;                  % uL
vjsr = 0.002*vcell;                 % uL
vcsr=0.008*vcell;                   % uL
vss  = 0.02*vcell;                  % uL
vsl=0.15*vcell;                     % uL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Reversal Potentials
ENa  = (R*T/F)*log(nao/nasl);
EK   = (R*T/F)*log(ko/ksl);
ECa= (R*T/(zca*F))*log(cao/casl);
PKNa = 0.01833;
EKs  = (R*T/F)*log((ko+PKNa*nao)/(ksl+PKNa*nasl));
%% Calcium_Fluxes_rate_constants from PRd
tautr1 = 120;
tautr2 = 120;	
gaptau = 12;
sstau = 0.2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CaMK Constants
KmCaMK = 0.15;  
aCaMK  = 0.05;  
bCaMK  = 0.00068;
CaMKo  = 0.05;  
KmCaM  = 0.0015;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% update CaMK -> X(41)
CaMKb = CaMKo*(1.0-CaMKt)/(1.0+KmCaM/cass);
CaMKa = CaMKb+CaMKt;
dCaMKt = aCaMK*CaMKb*(CaMKb+CaMKt) - bCaMK*CaMKt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                      MEMBRANE IONIC CURRENTS                        %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INa current from ORd
mss=1.0/ (1.0+exp((-(v+48.4264))/7.5653)); 
tm=1.0/(6.765*exp((v+11.64)/34.77)+8.552*exp(-(v+77.42)/5.955));
hss=1.0/(1+exp((v+78.5)/6.22)); 
thf=1.0/(3.6860e-06*exp(-(v+3.8875)/7.8579)+16*exp((v-0.4963)/9.1843));
ths=1.0/(0.009794*exp(-(v+17.95)/28.05)+0.3343*exp((v+5.730)/56.66));
Ahf=0.99;
Ahs=1.0-Ahf;
h=Ahf*hf+Ahs*hs;
jss=hss;
tj=4.8590+1.0/(0.8628*exp(-(v+116.7258)/7.6005)+1.1096*exp((v+6.2719)/9.0358)); 
hssp=1.0/(1+exp((v+84.7)/6.22));
thsp=3.0*ths;
hp=Ahf*hf+Ahs*hsp;
tjp=1.46*tj;
dm=(mss-m)/tm;
dhf=(hss-hf)/thf;
dhs=(hss-hs)/ths;
dj=(jss-j)/tj;
dhsp=(hssp-hsp)/thsp;
djp=(jss-jp)/tjp;
GNa=39.4572;
fINap=(1.0/(1.0+KmCaMK/CaMKa));
INa=GNa*(v-ENa)*m^3.0*((1.0-fINap)*h*j+fINap*hp*jp);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INaL current [from ORd]
mLss=1.0/(1.0+exp((-(v+42.85))/5.264));
tmL=tm;
hLss=1.0/(1.0+exp((v+87.61)/7.488));
thL=200.0;
hLssp=1.0/(1.0+exp((v+93.81)/7.488));
thLp=3.0*thL;
dmL=(mLss-mL)/tmL;
dhL=(hLss-hL)/thL;
dhLp=(hLssp-hLp)/thLp;
GNaL=0.0189; 
fINaLp=(1.0/(1.0+KmCaMK/CaMKa));
INaL=GNaL*(v-ENa)*mL*((1.0-fINaLp)*hL+fINaLp*hLp);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ICaL [from ORd]
vshift=15.19;
dss=1.0/(1.0+exp((-(v+3.940+3.3))/4.230));
td=0.6+1.0/(exp(-0.05*(v+6.0))+exp(0.09*(v+14.0)));
fss=1.0/(1.0+exp((v+19.58+3.3)/3.696));
tff=7.0+1.0/(0.0045*exp(-((v+vshift)+20.0)/10.0)+0.0045*exp(((v+vshift)+20.0)/10.0));
tfs=1000.0+1.0/(0.000035*exp(-((v+vshift)+5.0)/4.0)+0.000035*exp(((v+vshift)+5.0)/6.0)); 
Aff=0.6;
Afs=1.0-Aff;
f=Aff*ff+Afs*fs;
fcass=fss;
tfcaf=0.72*(7.0+1.0/(0.04*exp(-((v+vshift)-4.0)/7.0)+0.04*exp(((v+vshift)-4.0)/7.0)));
tfcas=0.49*(100.0 + 1.0/(0.00012*exp(-(v+vshift)/3.0)+0.00012*exp((v+vshift)/7.0)));
Afcaf=0.3+0.6/(1.0+exp((v-10.0)/10.0));
Afcas=1.0-Afcaf;
fca=Afcaf*fcaf+Afcas*fcas;
tjca=75.0;
ktaup=2.5;
tffp=ktaup*tff;
fp=Aff*ffp+Afs*fs;
tfcafp=ktaup*tfcaf;
fcap=Afcaf*fcafp+Afcas*fcas;
Kmn=0.002;
k2n=1000.0;
km2n=jca*1.0;
anca=1.0/(k2n/km2n+(1.0+Kmn/cass)^4.0);
dnca=(anca*k2n-nca*km2n);
PhiCaL=4.0*vffrt*(cass*exp(2.0*vfrt)-0.341*cao)/(exp(2.0*vfrt)-1.0);
PhiCaNa=1.0*vffrt*(0.75*nass*exp(1.0*vfrt)-0.75*nao)/(exp(1.0*vfrt)-1.0);
PhiCaK=1.0*vffrt*(0.75*kss*exp(1.0*vfrt)-0.75*ko)/(exp(1.0*vfrt)-1.0);
PCa=0.0001*0.776772893145948;
PCap=1.1*PCa;
PCaNa=0.00125*PCa;
PCaK=3.574e-4*PCa;
PCaNap=0.00125*PCap;
PCaKp=3.574e-4*PCap;
dd=(dss-d)/td;
dff=(fss-ff)/tff;
dfs=(fss-fs)/tfs;
dfcaf=(fcass-fcaf)/tfcaf;
dfcas=(fcass-fcas)/tfcas;
djca=(fcass-jca)/tjca;
dffp=(fss-ffp)/tffp;
dfcafp=(fcass-fcafp)/tfcafp;
fICaLp=(1.0/(1.0+KmCaMK/CaMKa));
ICaL=(1.0-fICaLp)*PCa*PhiCaL*d*(f*(1.0-nca)+jca*fca*nca)+...
             fICaLp*PCap*PhiCaL*d*(fp*(1.0-nca)+jca*fcap*nca);
ICaNa=(1.0-fICaLp)*PCaNa*PhiCaNa*d*(f*(1.0-nca)+jca*fca*nca)+...
                   fICaLp*PCaNap*PhiCaNa*d*(fp*(1.0-nca)+jca*fcap*nca);
ICaK=(1.0-fICaLp)*PCaK*PhiCaK*d*(f*(1.0-nca)+jca*fca*nca)+...
                  fICaLp*PCaKp*PhiCaK*d*(fp*(1.0-nca)+jca*fcap*nca);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% T_type_calcium_current_ICaT [from PRd]
gcat= 0.0754;      
bss= 1/(1+ exp (-(v+30)/7));
gss= 1/(1+exp((v+61)/5));
taub= 1/(1.068*exp((v+16.3)/30)+1.068*exp(-(v+16.3)/30));
taug= 1/(0.015*exp(-(v+71.7)/83.3)+0.015*exp((v+71.7)/15.4));
db= (bss-b)/taub;
dg= (gss-g)/taug;
ICaT= gcat*b*g*(v-ECa);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Ito current [from Han's data]
gtos= 0.192;          
ass= 1.0/(1.0+exp((20.0-v)/13.0));
iss= 1.0/(1.0+exp((v+27.0)/13.0));
atau= 1.0515/(1.0/(1.2089*(1.0+exp(-(v-18.4099)/29.3814)))+ 3.5/(1.0+exp((v+100.0)/29.3814))); 
tiS=43+1./(0.001416.*exp((-(v+96.52))/59.05)+1.780e-8.*exp((v+114.1)/8.079));
tiF=6.162+1./(0.3933.*exp((-(v+100.0))/100.0)+0.08004.*exp((v-8.0)/8.59));
da= (ass-a)/atau;
di= (iss-i)/tiS;
di2= (iss-i2)/tiF;
Ito= gtos*a*i*i2*(v-EK);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Sustained_outward_current [from Han's data]
g_sus=0.0301;
asus = 1.0/(1.0+exp(-(v-12)/16.0));
Isus= g_sus*asus*(v-EK);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% IKr current [modified from ORd]
xrss=1.0/(1.0+exp((-(v+8.337))/6.789));
txrf=12.98+1.0/(0.3652*exp((v-31.66+17.6)/3.869)+4.123e-5*exp((-(v-47.78+17.6))/20.38));
txrs=1.865+1.0/(0.06629*exp((v-34.70+17.2)/7.355)+1.128e-5*exp((-(v-29.74+17.2))/25.94));
Axrf=1.0/(1.0+exp((v+54.81)/38.21));
Axrs=1.0-Axrf;
dxrf=(xrss-xrf)/txrf;
dxrs=(xrss-xrs)/txrs;
xr=Axrf*xrf+Axrs*xrs;
rkr=1/(1+exp((v+55)/(0.32*75)))*1/(1+exp((v-10)/(0.32*30)));
GKr=0.0342;
IKr=GKr*sqrt(ko/5.4)*xr*rkr*(v-EK);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% IKs current [from ORd]
xs1ss=1.0/(1.0+exp((-(v+11.60))/8.932));
txs1=817.3+1.0/(2.326e-4*exp((v+48.28)/17.80)+0.001292*exp((-(v+210.0))/230.0));
dxs1=(xs1ss-xs1)/txs1;
xs2ss=xs1ss;
txs2=1.0/(0.01*exp((v-50.0)/20.0)+0.0193*exp((-(v+66.54))/31.0));
dxs2=(xs2ss-xs2)/txs2;
KsCa=1.0+0.6/(1.0+(3.8e-5/casl)^1.4);
GKs=0.0029;                       
IKs=GKs*KsCa*xs1*xs2*(v-EKs);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Hyper-polarization Activated Current (If) [From PRd]
yss= 1/(1+exp((v+87)/9.5));
tau_y= 2000/(exp(-(v+132)/10) + exp((v+57)/60));
dy=(yss-y)/tau_y;
IfNa  = 0.0116*y*y*(v-ENa);   
IfK = 0.0232*y*y*(v-EK);    
If = IfNa + IfK;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% IK1 current [from ORd's formulation with Han's data]
gk1=2.3238*sqrt(ko/5.4)*0.0455;               
xk1ss=1.0/(1.0+exp(-(v+2.5538*ko+144.59)/(1.5692*ko+3.8115)));
xk1tau=122.2/(exp((-(v+127.2))/20.36)+exp((v+236.8)/69.33));
dxk1=(xk1ss-xk1)/xk1tau;
rk1=1/(1.0+exp((v+116-5.5*ko)/11)); 
IK1=gk1*rk1*xk1*(v-EK);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INaCa current [from ORd]
kna1=15.0;      kna2=5.0;       kna3=88.12;     kasymm=12.5;
wna=6.0e4;      wca=6.0e4;      wnaca=5.0e3;    KmCaAct=150.0e-6;
kcaon=1.5e6;    kcaoff=5.0e3;   qna=0.5224;     qca=0.1670;
Gncx=9.5709e-04;        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INaCa_i current
hca=exp((qca*v*F)/(R*T));       hna=exp((qna*v*F)/(R*T));
h1=1+nasl/kna3*(1+hna);          h2=(nasl*hna)/(kna3*h1);
h3=1.0/h1;                      h4=1.0+nasl/kna1*(1+nasl/kna2);
h5=nasl*nasl/(h4*kna1*kna2);      h6=1.0/h4;
h7=1.0+nao/kna3*(1.0+1.0/hna);  h8=nao/(kna3*hna*h7);
h9=1.0/h7;                      h10=kasymm+1.0+nao/kna1*(1.0+nao/kna2);
h11=nao*nao/(h10*kna1*kna2);    h12=1.0/h10;

k1=h12*cao*kcaon;   k2=kcaoff;        k3p=h9*wca;     k3pp=h8*wnaca;      
k3=k3p+k3pp;        k4p=h3*wca/hca;   k4pp=h2*wnaca;  k4=k4p+k4pp;
k5=kcaoff;          k6=h6*casl*kcaon;  k7=h5*h2*wna;   k8=h8*h11*wna;

x1=k2*k4*(k7+k6)+k5*k7*(k2+k3); x2=k1*k7*(k4+k5)+k4*k6*(k1+k8);
x3=k1*k3*(k7+k6)+k8*k6*(k2+k3); x4=k2*k8*(k4+k5)+k3*k5*(k1+k8);

E1=x1/(x1+x2+x3+x4);    E2=x2/(x1+x2+x3+x4);
E3=x3/(x1+x2+x3+x4);    E4=x4/(x1+x2+x3+x4);

allo=1.0/(1.0+(KmCaAct/casl)^2.0);   
JncxNa=3.0*(E4*k7-E1*k8)+E3*k4pp-E2*k3pp;   
JncxCa=E2*k2-E1*k1;
%
INaCa_i=0.8*Gncx*allo*(zna*JncxNa+zca*JncxCa);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INaCa_ss current
h1=1+nass/kna3*(1+hna);         h2=(nass*hna)/(kna3*h1);
h3=1.0/h1;                      h4=1.0+nass/kna1*(1+nass/kna2);
h5=nass*nass/(h4*kna1*kna2);    h6=1.0/h4;
h7=1.0+nao/kna3*(1.0+1.0/hna);  h8=nao/(kna3*hna*h7);
h9=1.0/h7;                      h10=kasymm+1.0+nao/kna1*(1+nao/kna2);
h11=nao*nao/(h10*kna1*kna2);    h12=1.0/h10;

k1=h12*cao*kcaon;   k2=kcaoff;      k3p=h9*wca;     k3pp=h8*wnaca;
k3=k3p+k3pp;        k4p=h3*wca/hca; k4pp=h2*wnaca;  k4=k4p+k4pp;
k5=kcaoff;          k6=h6*cass*kcaon;   k7=h5*h2*wna;   k8=h8*h11*wna;

x1=k2*k4*(k7+k6)+k5*k7*(k2+k3);     x2=k1*k7*(k4+k5)+k4*k6*(k1+k8);
x3=k1*k3*(k7+k6)+k8*k6*(k2+k3);     x4=k2*k8*(k4+k5)+k3*k5*(k1+k8);

E1=x1/(x1+x2+x3+x4);    E2=x2/(x1+x2+x3+x4);
E3=x3/(x1+x2+x3+x4);    E4=x4/(x1+x2+x3+x4);

allo=1.0/(1.0+(KmCaAct/cass)^2.0);
JncxNa=3.0*(E4*k7-E1*k8)+E3*k4pp-E2*k3pp;
JncxCa=E2*k2-E1*k1;
%
INaCa_ss=0.2*Gncx*allo*(zna*JncxNa+zca*JncxCa);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INaK current [from ORd]
k1p=949.5;      k1m=182.4;      k2p=687.2;      k2m=39.4;
k3p=1899.0;     k3m=79300.0;    k4p=639.0;      k4m=40.0;
Knai0=9.073;    Knao0=27.78;    delta2=-0.1550;  
Knai=Knai0*exp((delta2*v*F)/(3.0*R*T));
Knao=Knao0*exp(((1.0-delta2)*v*F)/(3.0*R*T));
Kki=0.5;            Kko=0.3582;     MgADP=0.05;     MgATP=9.8;
Kmgatp=1.698e-7;    H=1.0e-7;       eP=4.2;         Khp=1.698e-7;   
Knap=224.0;         Kxkur=292.0;
P=eP/(1.0+H/Khp+nasl/Knap+ksl/Kxkur);

a1=(k1p*(nasl/Knai)^3.0)/((1.0+nasl/Knai)^3.0+(1.0+ksl/Kki)^2.0-1.0);
b1=k1m*MgADP;
a2=k2p;
b2=(k2m*(nao/Knao)^3.0)/((1.0+nao/Knao)^3.0+(1.0+ko/Kko)^2.0-1.0);
a3=(k3p*(ko/Kko)^2.0)/((1.0+nao/Knao)^3.0+(1.0+ko/Kko)^2.0-1.0);
b3=(k3m*P*H)/(1.0+MgATP/Kmgatp);
a4=(k4p*MgATP/Kmgatp)/(1.0+MgATP/Kmgatp);
b4=(k4m*(ksl/Kki)^2.0)/((1.0+nasl/Knai)^3.0+(1.0+ksl/Kki)^2.0-1.0);

x1=a4*a1*a2+b2*b4*b3+a2*b4*b3+b3*a1*a2;
x2=b2*b1*b4+a1*a2*a3+a3*b1*b4+a2*a3*b4;
x3=a2*a3*a4+b3*b2*b1+b2*b1*a4+a3*a4*b1;
x4=b4*b3*b2+a3*a4*a1+b2*a4*a1+b3*b2*a1;

E1=x1/(x1+x2+x3+x4);    E2=x2/(x1+x2+x3+x4);
E3=x3/(x1+x2+x3+x4);    E4=x4/(x1+x2+x3+x4);
JnakNa=3.0*(E1*a3-E2*b3);   JnakK=2.0*(E4*b1-E3*a1);    
Pnak=32.4872;
INaK=Pnak*(zna*JnakNa+zk*JnakK);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Background currents: INab, ICab [from ORd]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INab current
PNab = 3.75e-10*2.5;
INab = PNab*vffrt*(nasl*exp(vfrt)-nao)/(exp(vfrt)-1.0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ICab current
PCab = 2.5e-8;
ICab = PCab*4.0*vffrt*(casl*exp(2.0*vfrt)-0.341*cao)/(exp(2.0*vfrt)-1.0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% IpCa current  [from PRd]
ipcabar = 0.0005; 
kmpca = 0.0005;
IpCa= ipcabar/((kmpca/casl)+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Stimulation
amp = -40;      
duration = 1; %ms
if t <= duration
    Istim = amp;
else
    Istim = 0.0;
end
dv = - (INa+INaL+Ito+Isus+ICaL+ICaNa+ICaK+ICaT+If+IKr+IKs+IK1+INaCa_i+INaCa_ss+INaK+INab+IpCa+ICab+Istim);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                        SR_CA_FLUXES [from PRd]                      %%%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Diffusion Fluxes
JdiffNa= (nass-nasl)/sstau;
JgapNa=(nasl-nai)/gaptau;
JdiffK= (kss-ksl)/sstau;
JgapK=(ksl-ki)/gaptau;
Jdiff = (cass-casl)/sstau;  
Jgap  = (casl-cai)/gaptau;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Ip3R_Ca_Release
IP3 = 0.0001;  
k1 = 150000;
k1a = 16.5;
k0 = 96000;
k0a = 9.6;
k2 = 1800;
k2a = 0.21;
tauip3r = 3.7;
du=(cass*k2*(1-u)-k2a*u);
POip3 = tauip3r*IP3*cass*(1-u)/((1+IP3*k0/k0a)*(1+cass*k1/k1a));
Jip3 = 10.920*(cajsr-cass)*(POip3);                    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Ca_uptake_via_SERCA: qup1 & qup2
dqupcamkbar = 0.75;			        
dkmplbbar = 0.00017;	
kmup   = 0.00028;	
nsrbar = 15.0;	 
dkmplb		= dkmplbbar*CaMKa/(KmCaMK+CaMKa);
dqupcamk	= dqupcamkbar*CaMKa/(KmCaMK+CaMKa); 
Jup1		= (0.0002*(dqupcamk+1)/(1+((kmup-dkmplb)/casl))-0.00105*cansr/nsrbar);
Jup2		= (0.0026*(dqupcamk+1)/(1+((kmup-dkmplb)/cai))-0.0042*cansr/nsrbar);   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% RyR3_Ca_Release: qrel1 
REL1  = -((ICaL)*Acap/(vss*2.0*F)-(Jrel1 + Jip3)*vjsr/vss + Jdiff);
ireltau1 = 2*(1+1*(1/(1+((0.28/CaMKa)^8))))/(1+(0.0123/cajsr));
    if (REL1 > 0)
        irelss  = 15*(1+1*(1/(1+((0.28/CaMKa)^8))))*REL1/(1 + ((1.0/cajsr)^8));
    else
        irelss = 0;
    end
dJrel1= ((irelss-Jrel1)/ireltau1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% RyR2_Ca_release: qrel2 
REL2  = (-Jup2*vnsr/vmyo + Jgap*vsl/vmyo+ (Jrel2)*vcsr/vmyo);    
ireltau = 6*(1+1*(1/(1+((0.28/CaMKa)^8))))/(1+(0.0123/cacsr));
    if (REL2 > 0)
        irelss  = 91*(1+1*(1/(1+((0.28/CaMKa)^8))))*(REL2)/(1 + ((1/cacsr)^8)); 
    else
        irelss = 0;
    end
dJrel2= ((irelss-Jrel2)/ireltau);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Tranlocation Flux
Jtr1=(cansr-cajsr)/tautr1;
Jtr2=(cansr-cacsr)/tautr2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calcium Buffer Constants
%BSR
BSRmax=0.019975;
KmBSR=0.00087;
%BSL
BSLmax=0.4777;
KmBSL=0.0087;
%CSQN
csqnmax= 2.88;    
kmcsqn=0.8;     
csqnmax1 = 1.2;       
%CMDN
cmdnmax=0.1125;   
kmcmdn=0.00238;
cmdnmax1 = 1.25e-2;   
%TRPM
trpnmax=3.15e-2;
kmtrpn=0.0005;
trpnmax1 = 3.5e-3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                       IONIC CONCENTRATIONS                          %%%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% [Na]
dnass=-(ICaNa+3.0*INaCa_ss)*Acap/(F*vss*zna)-JdiffNa;
dnasl=-(INa+INaL+3.0*INaCa_i+3.0*INaK+IfNa+INab)*Acap/(F*vsl*zna)+JdiffNa*vss/vsl-JgapNa;
dnai=JgapNa*vsl/vmyo;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% [K]
dkss=-(ICaK)*Acap/(F*vss*zk)-JdiffK;
dksl=-(Ito+Isus+IKr+IKs+IfK+IK1+Istim-2.0*INaK)*Acap/(F*vsl*zk)+JdiffK*vss/vsl -JgapK;
dki=JgapK*vsl/vmyo;
%%%%%%%%%%%%%%%%%%%%%%%%% Intracellular [Ca] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% i (Myo)
trpnMYO=trpnmax*kmtrpn/(kmtrpn+cai)^2.0;
cmdnMYO=cmdnmax*kmcmdn/(kmcmdn+cai)^2.0;
Bcai   = 1.0 / (1.0+cmdnMYO+trpnMYO);
dcaibar= (-(Jup2)*vnsr/vmyo+Jgap*vsl/vmyo+(Jrel2)*vcsr/vmyo);
dcai = Bcai*dcaibar;
%% SS
Bcass  = 1.0/(1.0+BSRmax*KmBSR/(KmBSR+cass)^2.0+BSLmax*KmBSL/(KmBSL+cass)^2.0);
dcass  =  Bcass*(-(ICaL-2.0*INaCa_ss)*Acap/(2.0*F*vss)+(Jrel1+Jip3)*vjsr/vss-Jdiff);
%% SL
trpnSSL=(trpnmax1*kmtrpn)/((kmtrpn+casl)^2);
cmdnSSL=(cmdnmax1*kmcmdn)/((kmcmdn+casl)^2);
Bcasl=1/(1+trpnSSL+cmdnSSL);
dcaslbar= (-(Jup1)*vnsr/vsl+Jdiff*vss/vsl-Jgap-(IpCa+ICab+ICaT-2*INaCa_i)*Acap/(vsl*2.0*F));
dcasl=Bcasl*dcaslbar;
%%%%%%%%%%%%%%%%%%%%%%%%%% Sarcoplasmic [Ca] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NSR
dcansr= Jup1+Jup2-Jtr1*vjsr/vnsr-Jtr2*vcsr/vnsr;
%% JSR
csqn1=(csqnmax1*kmcsqn)/((kmcsqn+cajsr)^2);
Bcajsr = 1.0/(1.0+csqn1);
dcajsrbar=(Jtr1-Jrel1-Jip3);
dcajsr = Bcajsr*dcajsrbar;
%% CSR	
csqn=(csqnmax*kmcsqn)/((cacsr+kmcsqn)^2);
Bcacsr=1/(1+csqn);
dcacsrbar=(Jtr2-Jrel2);
dcacsr=Bcacsr*dcacsrbar;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                              OUTPUT                                 %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dX=[dv      dnai    dnasl   dnass   dki... 
    dksl    dkss    dcai    dcasl   dcass...
    dcansr  dcajsr  dcacsr  dm      dhf...
    dhs     dj      dhsp    djp     dmL...
    dhL     dhLp    da      di      di2...
    dd      dff     dfs     dfcaf   dfcas...
    djca    dnca    dffp    dfcafp  db...
    dg      dxrf    dxrs    dxs1    dxs2...
    dy      dxk1    dCaMKt  du      dJrel1...
    dJrel2]';  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end