- Author:
- Cristian Trovato <cristian.trovato@gmail.com>
- Date:
- 2020-04-07 19:26:00+01:00
- Desc:
- Model files upload
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/5d7/rawfile/95cadd4d554e8be258cf49a92eaf6c7fc4d66b96/Matlab/modTrovato2020.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
