print_lrd_web2007.m

Luo-Rudy cardiac cellular action potential model

function [currents]=print_LRd_web2007(State,t,Stm,data);


V=State(:,1);
H=State(:,2);
m=State(:,3);
J=State(:,4);
d=State(:,5);
f=State(:,6);
Xr=State(:,7);

Ca=conc_cai(State(:,8),data);
 JSR=conc_jsr(State(:,11),data);

Nai=State(:,9);
Ki=State(:,10);
%JSR=State(:,11);
NSR=State(:,12);
XS=State(:,13);
B=State(:,14);
G=State(:,15);
XS2=State(:,16);

%Qrelss=data.Qa./(1+(data.Qb./JSR).^8);

  Rel=State(:,17);
  Over=State(:,18);

%    fca=State(:,20);
   fca2 =1./(1+(Ca./data.kmca).^1);
currents.fca=fca2.*1	;
   EK = log(data.k_o./Ki)./data.frt; % potassium reverse potential
   ENa = log(data.na_o./Nai)./data.frt; % sodium reverse potential
   ECa = log(data.ca_o./Ca)./data.frt/2; % calcium reverse potential


   currents.ibarcal=(V*data.F*data.frt).*(data.gacai*Ca.*exp(2*V*data.frt)-data.gacao*data.ca_o)./(exp(2*V*data.frt)-1);
  currents.ical=d.*f.*data.pca*4.* currents.ibarcal.*currents.fca;
currents.icalna=d.*f.*data.pna.*(V*data.F*data.frt).*((data.ganai*Nai.*exp(V*data.frt)-data.ganao*data.na_o)./(exp(V*data.frt)-1)).*currents.fca;
 currents.icalk=d.*f.*data.pk.*(V*data.F*data.frt).*((data.gaki*Ki.*exp((V*data.frt))-data.gako*data.k_o)./(exp(V*data.frt)-1)).*currents.fca;
currents.icat =data.gcat.*B.*B.*G.*(V-ECa);

currents.INa=data.GNa*m.^3.*H.*J.*(V-ENa);
currents.inab = data.GNab*(V-ENa);
currents.G_Kr=data.gkrmax*sqrt(data.k_o/5.4)./(1+exp((V+9)/22.4)).*Xr;

currents.ikr =currents.G_Kr.*(V-EK);
currents.G_Ks= data.GKsmax*(1+data.IKsCa_max./(1+(data.IKsCa_Kd./Ca).^(1.4))).*XS.*XS2;
currents.iks = currents.G_Ks.*(V-log((data.k_o+data.prnak*data.na_o)./(Ki+data.prnak*Nai))/data.frt);
 currents.inaca = data.c1*exp((data.gammas-1)*V*data.frt).*(exp(V*data.frt).*Nai.^3*data.ca_o - data.na_o^3.*Ca)./...
     (1+data.c2*exp((data.gammas-1)*V*data.frt).*(exp(V*data.frt).*Nai.^3*data.ca_o + data.na_o^3*Ca));


currents.inak = data.ibarnak.*(1./(1+(data.kmnai./Nai).^2))*(1./(1+data.kmko./data.k_o))...
    ./(1+0.1245*exp((-0.1*V*data.frt))+0.0365*((exp(data.na_o/67.3)-1)/7).*exp(-V*data.frt));
currents.IKp = 0.00552.*(V-EK)./(1+exp((7.488-V)/5.98));


ak1 = 1.02./(1+exp(0.2385*(V-EK-59.215)));
bk1 = (0.49124*exp(0.08032*(V-EK+5.476))+exp(0.06175*(V-EK-594.31)))./...
    (1+exp(-0.5143*(V-EK+4.753)));
currents.G_K1= 0.75*sqrt(data.k_o/5.4)*ak1./(ak1+bk1);
currents.IK1 = currents.G_K1.*(V- EK);

%ikna = data.gkna.*0.85./(1+(data.kdkna./Nai).^(2.8)).*(0.8-0.65./(1+exp((V+125)./15))).*(V-EK);



currents.ipca = (data.ibarpca)./(1+data.kmpca./Ca);	 % sarcolema pump Ca SERCA

currents.ileak = data.iupbar./data.nsrbar*NSR;
%
currents.iup = data.iupbar./(1+(data.kmup./Ca).^1);
currents.itr = (NSR-JSR)./data.tautr;
currents.icab =data.gcab*(V- ECa);% background Ca
currents.csqn = data.csqnbar./(1+data.kmcsqn./JSR);
currents.trpn=data.trpnbar./(1+data.kmtrpn./Ca);
currents.cmdn=data.cmdnbar./(1+data.kmcmdn./Ca);
currents.naiont =    currents.INa+currents.inab+3*currents.inaca+currents.icalna+3*currents.inak;%insna+
currents.caiont = currents.ical+currents.icab+currents.ipca-2*currents.inaca+currents.icat;
currents.kiont = currents.ikr+currents.IK1+currents.IKp+currents.icalk-2*currents.inak+currents.iks-Stm';%+insk
%Over=0;
currents.dcai_out =-currents.caiont*data.AF/(data.vmyo*2);
currents.dcai_in=(currents.ileak-currents.iup)*data.vnsr/data.vmyo+(Over+Rel)*data.vjsr/data.vmyo;

currents.stm=Stm;

%currents.DVDT=  currents.naiont+ currents.caiont+currents.kiont;


currents.Cafree=Ca;
currents.JSRfree=JSR;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

















function [cai]=conc_cai(ca_t,data)
  	bmyo = data.cmdnbar+data.trpnbar-ca_t+data.kmtrpn+data.kmcmdn;
	cmyo = data.kmcmdn*data.kmtrpn -ca_t*(data.kmtrpn+data.kmcmdn)+data.trpnbar*data.kmcmdn+data.cmdnbar*data.kmtrpn;
	dmyo = -data.kmtrpn*data.kmcmdn*ca_t;
	


     cai = (2*(bmyo.*bmyo-3*cmyo).^(1/2)/3).*cos(acos((9*bmyo.*cmyo-2*bmyo.*bmyo.*bmyo-27*dmyo)./(2*(bmyo.*bmyo-3*cmyo).^1.5))/3)-(bmyo/3);


function [cajsr]=conc_jsr(ca_t,data)

b=data.csqnbar+data.kmcsqn-ca_t;
c=ca_t*data.kmcsqn;
cajsr=-b/2+(b.^2+4*c).^(1/2)/2;

function [cass]=conc_cass(ca_t,data)


	bmyo = data.BSLmax+data.BSRmax-ca_t+data.KmBSR+data.KmBSL;
	cmyo = data.KmBSL*data.KmBSR -ca_t*(data.KmBSR+data.KmBSL)+data.BSRmax*data.KmBSL+data.BSLmax*data.KmBSR;
	dmyo = -data.KmBSR*data.KmBSL*ca_t;
%      
     cass =( 2*(bmyo.*bmyo-3*cmyo).^(1/2)/3).*cos(acos((9*bmyo.*cmyo-2*bmyo.*bmyo.*bmyo-27*dmyo)./(2*(bmyo.*bmyo-3*cmyo).^1.5))/3)-(bmyo/3);
function [catrpn]=conc_catrpn(ca_t,trpn,data)
ca_t=ca_t-trpn;
b=data.cmdnbar+data.kmcmdn-ca_t;
c=ca_t*data.kmcmdn;
catrpn=-b/2+(b.^2+4*c).^(1/2)/2;