%     Cardiac model ToR-ORd-Land for acute and chronic myocardial infarction
%     Copyright (C) 2024 Xin Zhou. Contact: xin.zhou.joy@hotmail.com
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.
%
% The ToR-ORd model is described in the article "Development, calibration, and validation of a novel human ventricular myocyte model in health, disease, and drug block"
% by Jakub Tomek, Alfonso Bueno-Orovio, Elisa Passini, Xin Zhou, Ana Minchole, Oliver Britton, Chiara Bartolucci, Stefano Severi, Alvin Shrier, Laszlo Virag, Andras Varro and Blanca Rodriguez
% The article and supplemental materails are freely available in the
% Open Access jounal eLife
% Link to Article:
% https://elifesciences.org/articles/48890

% The models for acute and chronic infarction are described in the artile "Clinical phenotypes in acute and chronic infarction explained through human ventricular electromechanical modelling and simulations"
% by Xin Zhou, Zhinuo Jenny Wang, Julia Camps, Jakub Tomek, Alfonso Santiago, Adria Quintanas, Mariano Vazquez, Marmar Vaseghi and Blanca Rodriguez
% The article and supplemental materails are freely available in the
% Open Access jounal eLife
% Link to Article:
% https://elifesciences.org/reviewed-preprints/93002

% The model.m function contains the code for equations. The same models
% were implemented in other programming languages and simulation platforms,
% such as the the high-performance numerical software Alya and
% GPU-enabled open source MonoAlg3D simulator as shown in the following
% papers:
% https://elifesciences.org/reviewed-preprints/93002
% https://www.nature.com/articles/s41598-024-67951-5

function output=model(t,X,flag_ode,celltype,Parameters)

%%%%%%%% extract parameters from input
INa_Multiplier = Parameters(2);
INaL_Multiplier = Parameters(3);
Ito_Multiplier = Parameters(4);
ICaL_Multiplier = Parameters(5);
IKr_Multiplier = Parameters(6);
IKs_Multiplier = Parameters(7);
IK1_Multiplier = Parameters(8);
INaCa_Multiplier = Parameters(9);
INaK_Multiplier = Parameters(10);
Jrel_Multiplier = Parameters(11);
Jup_Multiplier = Parameters(12);
IKCa_Multiplier = Parameters(13);
ICaCl_Multiplier = Parameters(14);
CaMKPhosRate_Multiplier = Parameters(15);
Jrelptau_Multiplier = Parameters(16);
ICab_Multiplier = Parameters(17);

IKb_Multiplier = 1;
INab_Multiplier = 1;
IpCa_Multiplier = 1;
IClb_Multiplier = 1;

%%%%%% setup the ionic concentrations, ion channel distributions and the stimulus
%%%%%% current
nao = 140;
cao = 1.8;
ko = 5;
cli = 24;   % Intracellular Cl  [mM]
clo = 150;  % Extracellular Cl  [mM]
ICaL_fractionSS = 0.8;
INaCa_fractionSS = 0.35;
stimAmp = -53;
stimDur = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%give names to the state vector values
% voltage and ion concentrations
v=X(1);
nai=X(2);
nass=X(3);
ki=X(4);
kss=X(5);
cai=X(6);
cass=X(7);
cansr=X(8);
cajsr=X(9);
% INa
m=X(10);
hp=X(11);
h=X(12);
j=X(13);
jp=X(14);
% INaL
mL=X(15);
hL=X(16);
hLp=X(17);
% Ito
a=X(18);
iF=X(19);
iS=X(20);
ap=X(21);
iFp=X(22);
iSp=X(23);
% ICaL
d=X(24);
ff=X(25);
fs=X(26);
fcaf=X(27);
fcas=X(28);
jca=X(29);
nca=X(30);
nca_i=X(31);
ffp=X(32);
fcafp=X(33);
% IKs
xs1=X(34);
xs2=X(35);
% Jrel and CaMK
Jrelnp=X(36);
CaMKt=X(37);
% IKr
ikr_c0 = X(38);
ikr_c1 = X(39);
ikr_c2 = X(40);
ikr_o = X(41);
ikr_i = X(42);
% Jrel
Jrelp=X(43);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Land-Niederer model state variables for force geneation
XS=max(0,X(44));
XW=max(0,X(45));
Ca_TRPN=max(0,X(46));
TmBlocked=X(47);
ZETAS=X(48);
ZETAW=X(49);
%% Land-Niederer model
%==========================================================================
% input coded here:
mode='intact';
lambda=1;
lambda_rate=0;


%EC parameters
perm50=0.35;
TRPN_n=2;
koff=0.1;
dr=0.25;
wfrac=0.5;
TOT_A=25;
ktm_unblock= 0.021;
beta_1=-2.4;
beta_0=2.3;
gamma=0.0085;
gamma_wu=0.615;
phi=2.23;

if strcmp(mode,'skinned')
  nperm=2.2;
  ca50=2.5;
  Tref=40.5;
  nu=1;
  mu=1;
else
  nperm=2.036;
  ca50=0.805;
  Tref=120;
  nu=7;
  mu=3;
end

k_ws=0.004;
k_uw=0.026;

lambda_min=0.87;
lambda_max=1.2;

dydt=zeros(1,6);

k_ws=k_ws*mu;
k_uw=k_uw*nu;

cdw=phi*k_uw*(1-dr)*(1-wfrac)/((1-dr)*wfrac);
cds=phi*k_ws*(1-dr)*wfrac/dr;
k_wu=k_uw*(1/wfrac-1)-k_ws;
k_su=k_ws*(1/dr-1)*wfrac;
A=(0.25*TOT_A)/((1-dr)*wfrac+dr)*(dr/0.25);

%XB model
lambda0=min(lambda_max,lambda);
Lfac=max(0,1+beta_0*(lambda0+min(lambda_min,lambda0)-(1+lambda_min)));

XU=(1-TmBlocked)-XW-XS; % unattached available xb = all - tm blocked - already prepowerstroke - already post-poststroke - no overlap
xb_ws=k_ws*XW;
xb_uw=k_uw*XU ;
xb_wu=k_wu*XW;
xb_su=k_su*XS;

gamma_rate=gamma*max((ZETAS>0).*ZETAS,(ZETAS<-1).*(-ZETAS-1));
xb_su_gamma=gamma_rate*XS;
gamma_rate_w=gamma_wu*abs(ZETAW); % weak xbs don't like being strained
xb_wu_gamma=gamma_rate_w*XW;

dydt(1)=xb_ws-xb_su-xb_su_gamma;
dydt(2)=xb_uw-xb_wu-xb_ws-xb_wu_gamma;

ca50=ca50+beta_1*min(0.2,lambda-1);
dydt(3)=koff*(((cai*1000)/ca50)^TRPN_n*(1-Ca_TRPN)-Ca_TRPN); % untouched

XSSS=dr*0.5;
XWSS=(1-dr)*wfrac*0.5;
ktm_block=ktm_unblock*(perm50^nperm)*0.5/(0.5-XSSS-XWSS);
dydt(4)=ktm_block*min(100,(Ca_TRPN^-(nperm/2)))*XU-ktm_unblock*(Ca_TRPN^(nperm/2))*TmBlocked;

%velocity dependence -- assumes distortion resets on W->S
dydt(5)=A*lambda_rate-cds*ZETAS;% - gamma_rate * ZETAS;
dydt(6)=A*lambda_rate-cdw*ZETAW;% - gamma_rate_w * ZETAW;

% Active Force
Ta=Lfac*(Tref/dr)*((ZETAS+1)*XS+(ZETAW)*XW);
%========================================================================


%cell geometry
L=0.01;
rad=0.0011;
vcell=1000*3.14*rad*rad*L;
Ageo=2*3.14*rad*rad+2*3.14*rad*L;
Acap=2*Ageo;
vmyo=0.68*vcell;
vnsr=0.0552*vcell;
vjsr=0.0048*vcell;
vss=0.02*vcell;

%CaMK constants
KmCaMK=0.15;

aCaMK=0.05*CaMKPhosRate_Multiplier;
bCaMK=0.00068;
CaMKo=0.05;
KmCaM=0.0015;
%update CaMK
CaMKb=CaMKo*(1.0-CaMKt)/(1.0+KmCaM/cass);
CaMKa=CaMKb+CaMKt;
dCaMKt=aCaMK*CaMKb*(CaMKb+CaMKt)-bCaMK*CaMKt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%physical constants
R=8314.0;
T=310.0;
F=96485.0;

%reversal potentials
ENa=(R*T/F)*log(nao/nai);
EK=(R*T/F)*log(ko/ki);
PKNa=0.01833;
EKs=(R*T/F)*log((ko+PKNa*nao)/(ki+PKNa*nai));

%convenient shorthand calculations
vffrt=v*F*F/(R*T);
vfrt=v*F/(R*T);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% fractions of ion channel phosphorylation
fINap=(1.0/(1.0+KmCaMK/CaMKa));
fINaLp=(1.0/(1.0+KmCaMK/CaMKa));
fItop=(1.0/(1.0+KmCaMK/CaMKa));
fICaLp=(1.0/(1.0+KmCaMK/CaMKa));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INa formulations
% The Grandi implementation updated with INa phosphorylation.
% m gate
mss = 1 / ((1 + exp( -(56.86 + v) / 9.03 ))^2);
taum = 0.1292 * exp(-((v+45.79)/15.54)^2) + 0.06487 * exp(-((v-4.823)/51.12)^2);
dm = (mss - m) / taum;

% h gate
ah = (v >= -40) * (0)...
    + (v < -40) * (0.057 * exp( -(v + 80) / 6.8 ));
bh = (v >= -40) * (0.77 / (0.13*(1 + exp( -(v + 10.66) / 11.1 )))) ...
    + (v < -40) * ((2.7 * exp( 0.079 * v) + 3.1*10^5 * exp(0.3485 * v)));
tauh = 1 / (ah + bh);
hss = 1 / ((1 + exp( (v + 71.55)/7.43 ))^2);
dh = (hss - h) / tauh;
% j gate
aj = (v >= -40) * (0) ...
    +(v < -40) * (((-2.5428 * 10^4*exp(0.2444*v) - 6.948*10^-6 * exp(-0.04391*v)) * (v + 37.78)) / ...
    (1 + exp( 0.311 * (v + 79.23) )));
bj = (v >= -40) * ((0.6 * exp( 0.057 * v)) / (1 + exp( -0.1 * (v + 32) ))) ...
    + (v < -40) * ((0.02424 * exp( -0.01052 * v )) / (1 + exp( -0.1378 * (v + 40.14) )));
tauj = 1 / (aj + bj);
jss = 1 / ((1 + exp( (v + 71.55)/7.43 ))^2);
dj = (jss - j) / tauj;

% h phosphorylated
hssp = 1 / ((1 + exp( (v + 71.55 + 6)/7.43 ))^2);
dhp = (hssp - hp) / tauh;
% j phosphorylated
taujp = 1.46 * tauj;
djp = (jss - jp) / taujp;

GNa = 11.7802;
INa=INa_Multiplier * GNa*(v-ENa)*m^3.0*((1.0-fINap)*h*j+fINap*hp*jp);

%% INaL formulation
%calculate INaL
mLss=1.0/(1.0+exp((-(v+42.85))/5.264));
tm = 0.1292 * exp(-((v+45.79)/15.54)^2) + 0.06487 * exp(-((v-4.823)/51.12)^2);
tmL=tm;
dmL=(mLss-mL)/tmL;
hLss=1.0/(1.0+exp((v+87.61)/7.488));
thL=200.0;
dhL=(hLss-hL)/thL;
hLssp=1.0/(1.0+exp((v+93.81)/7.488));
thLp=3.0*thL;
dhLp=(hLssp-hLp)/thLp;
GNaL=0.0279 * INaL_Multiplier;
if celltype==1
    GNaL=GNaL*0.6;
end

INaL=GNaL*(v-ENa)*mL*((1.0-fINaLp)*hL+fINaLp*hLp);


%% Ito formulation
%calculate Ito
ass=1.0/(1.0+exp((-(v-14.34))/14.82));
ta=1.0515/(1.0/(1.2089*(1.0+exp(-(v-18.4099)/29.3814)))+3.5/(1.0+exp((v+100.0)/29.3814)));
da=(ass-a)/ta;
iss=1.0/(1.0+exp((v+43.94)/5.711));
if celltype==1
    delta_epi=1.0-(0.95/(1.0+exp((v+70.0)/5.0)));
else
    delta_epi=1.0;
end
tiF=4.562+1/(0.3933*exp((-(v+100.0))/100.0)+0.08004*exp((v+50.0)/16.59));
tiS=23.62+1/(0.001416*exp((-(v+96.52))/59.05)+1.780e-8*exp((v+114.1)/8.079));
tiF=tiF*delta_epi;
tiS=tiS*delta_epi;
AiF=1.0/(1.0+exp((v-213.6)/151.2));
AiS=1.0-AiF;
diF=(iss-iF)/tiF;
diS=(iss-iS)/tiS;
i=AiF*iF+AiS*iS;
assp=1.0/(1.0+exp((-(v-24.34))/14.82));
dap=(assp-ap)/ta;
dti_develop=1.354+1.0e-4/(exp((v-167.4)/15.89)+exp(-(v-12.23)/0.2154));
dti_recover=1.0-0.5/(1.0+exp((v+70.0)/20.0));
tiFp=dti_develop*dti_recover*tiF;
tiSp=dti_develop*dti_recover*tiS;
diFp=(iss-iFp)/tiFp;
diSp=(iss-iSp)/tiSp;
ip=AiF*iFp+AiS*iSp;
Gto=0.16 * Ito_Multiplier;
if celltype==1
    Gto=Gto*2.0;
elseif celltype==2
    Gto=Gto*2.0;
end

Ito=Gto*(v-EK)*((1.0-fItop)*a*i+fItop*ap*ip);


%% ICaL formulation
%calculate ICaL, ICaNa, ICaK

dss=1.0763*exp(-1.0070*exp(-0.0829*(v)));  % magyar
if(v >31.4978) % activation cannot be greater than 1
    dss = 1;
end


td= 0.6+1.0/(exp(-0.05*(v+6.0))+exp(0.09*(v+14.0)));

dd=(dss-d)/td;
fss=1.0/(1.0+exp((v+19.58)/3.696));
tff=7.0+1.0/(0.0045*exp(-(v+20.0)/10.0)+0.0045*exp((v+20.0)/10.0));
tfs=1000.0+1.0/(0.000035*exp(-(v+5.0)/4.0)+0.000035*exp((v+5.0)/6.0));
Aff=0.6;
Afs=1.0-Aff;
dff=(fss-ff)/tff;
dfs=(fss-fs)/tfs;
f=Aff*ff+Afs*fs;
fcass=fss;
tfcaf=7.0+1.0/(0.04*exp(-(v-4.0)/7.0)+0.04*exp((v-4.0)/7.0));
tfcas=100.0+1.0/(0.00012*exp(-v/3.0)+0.00012*exp(v/7.0));

Afcaf=0.3+0.6/(1.0+exp((v-10.0)/10.0));

Afcas=1.0-Afcaf;
dfcaf=(fcass-fcaf)/tfcaf;
dfcas=(fcass-fcas)/tfcas;
fca=Afcaf*fcaf+Afcas*fcas;

tjca = 75;
jcass = 1.0/(1.0+exp((v+18.08)/(2.7916)));
djca=(jcass-jca)/tjca;
tffp=2.5*tff;
dffp=(fss-ffp)/tffp;
fp=Aff*ffp+Afs*fs;
tfcafp=2.5*tfcaf;
dfcafp=(fcass-fcafp)/tfcafp;
fcap=Afcaf*fcafp+Afcas*fcas;

% SS nca
Kmn=0.002;
k2n=500.0;
km2n=jca*1;
anca=1.0/(k2n/km2n+(1.0+Kmn/cass)^4.0);
dnca=anca*k2n-nca*km2n;

% myoplasmic nca
anca_i = 1.0/(k2n/km2n+(1.0+Kmn/cai)^4.0);
dnca_i = anca_i*k2n-nca_i*km2n;

% SS driving force

Io = 0.5*(nao + ko + clo + 4*cao)/1000 ; % ionic strength outside. /1000 is for things being in micromolar
Ii = 0.5*(nass + kss + cli + 4*cass)/1000 ; % ionic strength outside. /1000 is for things being in micromolar
% The ionic strength is too high for basic DebHuc. We'll use Davies
dielConstant = 74; % water at 37�.
temp = 310; % body temp in kelvins.
constA = 1.82*10^6*(dielConstant*temp)^(-1.5);

gamma_cai = exp(-constA * 4 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii));
gamma_cao = exp(-constA * 4 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io));
gamma_nai = exp(-constA * 1 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii));
gamma_nao = exp(-constA * 1 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io));
gamma_ki = exp(-constA * 1 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii));
gamma_kao = exp(-constA * 1 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io));


PhiCaL_ss =  4.0*vffrt*(gamma_cai*cass*exp(2.0*vfrt)-gamma_cao*cao)/(exp(2.0*vfrt)-1.0);
PhiCaNa_ss =  1.0*vffrt*(gamma_nai*nass*exp(1.0*vfrt)-gamma_nao*nao)/(exp(1.0*vfrt)-1.0);
PhiCaK_ss =  1.0*vffrt*(gamma_ki*kss*exp(1.0*vfrt)-gamma_kao*ko)/(exp(1.0*vfrt)-1.0);

% Myo driving force
Io = 0.5*(nao + ko + clo + 4*cao)/1000 ; % ionic strength outside. /1000 is for things being in micromolar
Ii = 0.5*(nai + ki + cli + 4*cai)/1000 ; % ionic strength outside. /1000 is for things being in micromolar
% The ionic strength is too high for basic DebHuc. We'll use Davies
dielConstant = 74; % water at 37�.
temp = 310; % body temp in kelvins.
constA = 1.82*10^6*(dielConstant*temp)^(-1.5);

gamma_cai = exp(-constA * 4 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii));
gamma_cao = exp(-constA * 4 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io));
gamma_nai = exp(-constA * 1 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii));
gamma_nao = exp(-constA * 1 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io));
gamma_ki = exp(-constA * 1 * (sqrt(Ii)/(1+sqrt(Ii))-0.3*Ii));
gamma_kao = exp(-constA * 1 * (sqrt(Io)/(1+sqrt(Io))-0.3*Io));

gammaCaoMyo = gamma_cao;
gammaCaiMyo = gamma_cai;

PhiCaL_i =  4.0*vffrt*(gamma_cai*cai*exp(2.0*vfrt)-gamma_cao*cao)/(exp(2.0*vfrt)-1.0);
PhiCaNa_i =  1.0*vffrt*(gamma_nai*nai*exp(1.0*vfrt)-gamma_nao*nao)/(exp(1.0*vfrt)-1.0);
PhiCaK_i =  1.0*vffrt*(gamma_ki*ki*exp(1.0*vfrt)-gamma_kao*ko)/(exp(1.0*vfrt)-1.0);
% The rest
PCa=8.3757e-05 * ICaL_Multiplier;

if celltype==1
    PCa=PCa*1.2;
elseif celltype==2
    PCa=PCa*1.8; %2;
end

PCap=1.1*PCa;
PCaNa=0.00125*PCa;
PCaK=3.574e-4*PCa;
PCaNap=0.00125*PCap;
PCaKp=3.574e-4*PCap;

ICaL_ss=(1.0-fICaLp)*PCa*PhiCaL_ss*d*(f*(1.0-nca)+jca*fca*nca)+fICaLp*PCap*PhiCaL_ss*d*(fp*(1.0-nca)+jca*fcap*nca);
ICaNa_ss=(1.0-fICaLp)*PCaNa*PhiCaNa_ss*d*(f*(1.0-nca)+jca*fca*nca)+fICaLp*PCaNap*PhiCaNa_ss*d*(fp*(1.0-nca)+jca*fcap*nca);
ICaK_ss=(1.0-fICaLp)*PCaK*PhiCaK_ss*d*(f*(1.0-nca)+jca*fca*nca)+fICaLp*PCaKp*PhiCaK_ss*d*(fp*(1.0-nca)+jca*fcap*nca);

ICaL_i=(1.0-fICaLp)*PCa*PhiCaL_i*d*(f*(1.0-nca_i)+jca*fca*nca_i)+fICaLp*PCap*PhiCaL_i*d*(fp*(1.0-nca_i)+jca*fcap*nca_i);
ICaNa_i=(1.0-fICaLp)*PCaNa*PhiCaNa_i*d*(f*(1.0-nca_i)+jca*fca*nca_i)+fICaLp*PCaNap*PhiCaNa_i*d*(fp*(1.0-nca_i)+jca*fcap*nca_i);
ICaK_i=(1.0-fICaLp)*PCaK*PhiCaK_i*d*(f*(1.0-nca_i)+jca*fca*nca_i)+fICaLp*PCaKp*PhiCaK_i*d*(fp*(1.0-nca_i)+jca*fcap*nca_i);


% And we weight ICaL (in ss) and ICaL_i
ICaL_i = ICaL_i * (1-ICaL_fractionSS);
ICaNa_i = ICaNa_i * (1-ICaL_fractionSS);
ICaK_i = ICaK_i * (1-ICaL_fractionSS);
ICaL_ss = ICaL_ss * ICaL_fractionSS;
ICaNa_ss = ICaNa_ss * ICaL_fractionSS;
ICaK_ss = ICaK_ss * ICaL_fractionSS;

ICaL = ICaL_ss + ICaL_i;
ICaNa = ICaNa_ss + ICaNa_i;
ICaK = ICaK_ss + ICaK_i;
ICaL_tot = ICaL + ICaNa + ICaK;

%% IKr formulation
%calculate IKr
% transition rates
% from c0 to c1 in l-v model,
alpha = 0.1161 * exp(0.2990 * vfrt);
% from c1 to c0 in l-v/
beta =  0.2442 * exp(-1.604 * vfrt);

% from c1 to c2 in l-v/
alpha1 = 1.25 * 0.1235 ;
% from c2 to c1 in l-v/
beta1 =  0.1911;

% from c2 to o/           c1 to o
alpha2 =0.0578 * exp(0.9710 * vfrt); %
% from o to c2/
beta2 = 0.349e-3* exp(-1.062 * vfrt); %

% from o to i
alphai = 0.2533 * exp(0.5953 * vfrt); %
% from i to o
betai = 1.25* 0.0522 * exp(-0.8209 * vfrt); %

% from c2 to i (from c1 in orig)
alphac2ToI = 0.52e-4 * exp(1.525 * vfrt); %
% from i to c2
% betaItoC2 = 0.85e-8 * exp(-1.842 * vfrt); %
betaItoC2 = (beta2 * betai * alphac2ToI)/(alpha2 * alphai); %
% transitions themselves
% for reason of backward compatibility of naming of an older version of a
% MM IKr, c3 in code is c0 in article diagram, c2 is c1, c1 is c2.

dt_ikr_c0 = ikr_c1 * beta - ikr_c0 * alpha; % delta for c0
dt_ikr_c1 = ikr_c0 * alpha + ikr_c2*beta1 - ikr_c1*(beta+alpha1); % c1
dt_ikr_c2 = ikr_c1 * alpha1 + ikr_o*beta2 + ikr_i*betaItoC2 - ikr_c2 * (beta1 + alpha2 + alphac2ToI); % subtraction is into c2, to o, to i. % c2
dt_ikr_o = ikr_c2 * alpha2 + ikr_i*betai - ikr_o*(beta2+alphai);
dt_ikr_i = ikr_c2*alphac2ToI + ikr_o*alphai - ikr_i*(betaItoC2 + betai);

GKr = 0.0321 * sqrt(ko/5) * IKr_Multiplier; % 1st element compensates for change to ko (sqrt(5/5.4)* 0.0362)
if celltype==1
    GKr=GKr*1.3;
elseif celltype==2
    GKr=GKr*0.8;
end

IKr = GKr * ikr_o  * (v-EK);

%% IKs formulation
%calculate IKs
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/cai)^1.4);
GKs= 0.0011 * IKs_Multiplier;
if celltype==1
    GKs=GKs*1.4;
end
IKs=GKs*KsCa*xs1*xs2*(v-EKs);

%% IK1 formulation
% IK1
aK1 = 4.094/(1+exp(0.1217*(v-EK-49.934)));
bK1 = (15.72*exp(0.0674*(v-EK-3.257))+exp(0.0618*(v-EK-594.31)))/(1+exp(-0.1629*(v-EK+14.207)));
K1ss = aK1/(aK1+bK1);

GK1=IK1_Multiplier  * 0.6992; %0.7266; %* sqrt(5/5.4))
if celltype==1
    GK1=GK1*1.2;
elseif celltype==2
    GK1=GK1*1.3;
end
IK1=GK1*sqrt(ko/5)*K1ss*(v-EK);

%% INaCa
%calculate INaCa_i
zca = 2.0;
kna1=15.0;
kna2=5.0;
kna3=88.12;
kasymm=12.5;
wna=6.0e4;
wca=6.0e4;
wnaca=5.0e3;
kcaon=1.5e6;
kcaoff=5.0e3;
qna=0.5224;
qca=0.1670;
hca=exp((qca*v*F)/(R*T));
hna=exp((qna*v*F)/(R*T));
h1=1+nai/kna3*(1+hna);
h2=(nai*hna)/(kna3*h1);
h3=1.0/h1;
h4=1.0+nai/kna1*(1+nai/kna2);
h5=nai*nai/(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*cai*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);
KmCaAct=150.0e-6;
allo=1.0/(1.0+(KmCaAct/cai)^2.0);
zna=1.0;
JncxNa=3.0*(E4*k7-E1*k8)+E3*k4pp-E2*k3pp;
JncxCa=E2*k2-E1*k1;
Gncx= 0.0034* INaCa_Multiplier;
if celltype==1
    Gncx=Gncx*1.1;
elseif celltype==2
    Gncx=Gncx*1.4;
end
INaCa_i=(1-INaCa_fractionSS)*Gncx*allo*(zna*JncxNa+zca*JncxCa);

%calculate INaCa_ss
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);
KmCaAct=150.0e-6 ;
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=INaCa_fractionSS*Gncx*allo*(zna*JncxNa+zca*JncxCa);

%% INaK formulation
%calculate INaK
zna=1.0;
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;
delta=-0.1550;
Knai=Knai0*exp((delta*v*F)/(3.0*R*T));
Knao=Knao0*exp(((1.0-delta)*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+nai/Knap+ki/Kxkur);
a1=(k1p*(nai/Knai)^3.0)/((1.0+nai/Knai)^3.0+(1.0+ki/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*(ki/Kki)^2.0)/((1.0+nai/Knai)^3.0+(1.0+ki/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);
zk=1.0;
JnakNa=3.0*(E1*a3-E2*b3);
JnakK=2.0*(E4*b1-E3*a1);
Pnak= 15.4509 * INaK_Multiplier;
if celltype==1
    Pnak=Pnak*0.9;
elseif celltype==2
    Pnak=Pnak*0.7;
end

INaK=Pnak*(zna*JnakNa+zk*JnakK);

%% IKb formulation
%calculate IKb
xkb=1.0/(1.0+exp(-(v-10.8968)/(23.9871)));
GKb=0.0189*IKb_Multiplier*0.9;
if celltype==1
    GKb=GKb*0.6;
end
IKb=GKb*xkb*(v-EK);

%% INab formulation
%calculate INab
PNab=1.9239e-09*INab_Multiplier;
INab=PNab*vffrt*(nai*exp(vfrt)-nao)/(exp(vfrt)-1.0);

%% ICab formulation
%calculate ICab
PCab=5.9194e-08*ICab_Multiplier;
ICab=PCab*4.0*vffrt*(gammaCaiMyo*cai*exp(2.0*vfrt)-gammaCaoMyo*cao)/(exp(2.0*vfrt)-1.0);

%% IpCa formulation
%calculate IpCa
GpCa=5e-04*IpCa_Multiplier;
IpCa=GpCa*cai/(0.0005+cai);

%% IKCa formulation
% calcuate Ca activated K current, IKCa (SK)
gkca= 0.003;
ikcan=3.5;
kdikca=6.05e-04;
FractionIKCass=0.8;
FractionIKCai=1-FractionIKCass;
IKCa_ss=gkca*IKCa_Multiplier*FractionIKCass*(cass^ikcan)/(cass^ikcan+kdikca^ikcan)*(v-EK);
IKCa_i=gkca*IKCa_Multiplier*FractionIKCai*(cai^ikcan)/(cai^ikcan+kdikca^ikcan)*(v-EK);
IKCa=IKCa_ss+IKCa_i;

%% IClCa formulation
% I_ClCa: Ca-activated Cl Current, I_Clbk: background Cl Current

ecl = (R*T/F)*log(cli/clo);            % [mV]

Fjunc = 1;   Fsl = 1-Fjunc; % fraction in SS and in myoplasm - as per literature, I(Ca)Cl is in junctional subspace


GClCa = ICaCl_Multiplier * 0.2843;   % [mS/uF]
GClB = IClb_Multiplier * 1.98e-3;        % [mS/uF] %
KdClCa = 0.1;    % [mM]

I_ClCa_junc = Fjunc*GClCa/(1+KdClCa/cass)*(v-ecl);
I_ClCa_sl = Fsl*GClCa/(1+KdClCa/cai)*(v-ecl);

I_ClCa = I_ClCa_junc+I_ClCa_sl;
I_Clbk = GClB*(v-ecl);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% Jrel formulation
%calculate ryanodione receptor calcium induced calcium release from the jsr
fJrelp=(1.0/(1.0+KmCaMK/CaMKa));

jsrMidpoint = 1.7;

bt=4.75;
a_rel=0.5*bt;
Jrel_inf=a_rel*(-ICaL_ss)/(1.0+(jsrMidpoint/cajsr)^8.0);
if celltype==2
    Jrel_inf=Jrel_inf*1.7;
end
tau_rel=bt/(1.0+0.0123/cajsr);

if tau_rel<0.001
    tau_rel=0.001;
end

dJrelnp=(Jrel_inf-Jrelnp)/tau_rel;
btp=1.25*bt;
a_relp=0.5*btp;
Jrel_infp=a_relp*(-ICaL_ss)/(1.0+(jsrMidpoint/cajsr)^8.0);
if celltype==2
    Jrel_infp=Jrel_infp*1.7;
end
tau_relp=Jrelptau_Multiplier*btp/(1.0+0.0123/cajsr);

if tau_relp<0.001
    tau_relp=0.001;
end

dJrelp=(Jrel_infp-Jrelp)/tau_relp;

Jrel=Jrel_Multiplier * 1.5378 * ((1.0-fJrelp)*Jrelnp+fJrelp*Jrelp);

%% Jup formulation
fJupp=(1.0/(1.0+KmCaMK/CaMKa));
%calculate serca pump, ca uptake flux
Jupnp=Jup_Multiplier * 0.005425*cai/(cai+0.00092);
Jupp=Jup_Multiplier * 2.75*0.005425*cai/(cai+0.00092-0.00017);
if celltype==1
    Jupnp=Jupnp*1.3;
    Jupp=Jupp*1.3;
end

Jleak=Jup_Multiplier* 0.0048825*cansr/15.0;
Jup=(1.0-fJupp)*Jupnp+fJupp*Jupp-Jleak;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Jtr=(cansr-cajsr)/60;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% calculate the stimulus current, Istim
amp=stimAmp;
duration=stimDur;
if t<=duration
    Istim=amp;
else
    Istim=0.0;
end
%%
%update the membrane voltage

dv=-(INa+INaL+Ito+ICaL+ICaNa+ICaK+IKr+IKs+IK1+INaCa_i+INaCa_ss+INaK+INab+IKb+IpCa+ICab+I_ClCa+I_Clbk+Istim+IKCa);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%calculate diffusion fluxes
JdiffNa=(nass-nai)/2.0;
JdiffK=(kss-ki)/2.0;
Jdiff=(cass-cai)/0.2;

%calcium buffer constants
cmdnmax= 0.05;
if celltype==1
    cmdnmax=cmdnmax*1.3;
end
kmcmdn=0.00238;
trpnmax=0.07;
% kmtrpn=0.0005;
BSRmax=0.047;
KmBSR = 0.00087;
BSLmax=1.124;
KmBSL = 0.0087;
csqnmax=10.0;
kmcsqn=0.8;

%update intracellular concentrations, using buffers for cai, cass, cajsr
dnai=-(ICaNa_i+INa+INaL+3.0*INaCa_i+3.0*INaK+INab)*Acap/(F*vmyo)+JdiffNa*vss/vmyo;
dnass=-(ICaNa_ss+3.0*INaCa_ss)*Acap/(F*vss)-JdiffNa;

dki=-(ICaK_i+Ito+IKr+IKs+IK1+IKb+Istim-2.0*INaK+IKCa_i)*Acap/(F*vmyo)+JdiffK*vss/vmyo;
dkss=-(ICaK_ss+IKCa_ss)*Acap/(F*vss)-JdiffK;


Bcai=1.0/(1.0+cmdnmax*kmcmdn/(kmcmdn+cai)^2.0);
dcai=Bcai*(-(ICaL_i + IpCa+ICab-2.0*INaCa_i)*Acap/(2.0*F*vmyo)-Jup*vnsr/vmyo+Jdiff*vss/vmyo-dydt(3)*trpnmax);

Bcass=1.0/(1.0+BSRmax*KmBSR/(KmBSR+cass)^2.0+BSLmax*KmBSL/(KmBSL+cass)^2.0);
dcass=Bcass*(-(ICaL_ss-2.0*INaCa_ss)*Acap/(2.0*F*vss)+Jrel*vjsr/vss-Jdiff);

dcansr=Jup-Jtr*vjsr/vnsr;

Bcajsr=1.0/(1.0+csqnmax*kmcsqn/(kmcsqn+cajsr)^2.0);
dcajsr=Bcajsr*(Jtr-Jrel);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%output the state venctor when ode_flag==1, and the calculated currents and fluxes otherwise
if flag_ode==1
    output=[dv dnai dnass dki dkss dcai dcass dcansr dcajsr dm dhp dh dj djp dmL dhL dhLp da diF diS dap diFp diSp,...
        dd dff dfs dfcaf dfcas djca dnca dnca_i dffp dfcafp dxs1 dxs2 dJrelnp dCaMKt,...
        dt_ikr_c0 dt_ikr_c1 dt_ikr_c2 dt_ikr_o dt_ikr_i dJrelp dydt(1) dydt(2) dydt(3) dydt(4) dydt(5) dydt(6)]';
else
    output=[INa INaL Ito ICaL IKr IKs IK1 INaCa_i INaCa_ss INaK  IKb INab ICab IpCa Jdiff JdiffNa JdiffK Jup Jleak Jtr Jrel CaMKa Istim, fINap, ...
        fINaLp, fICaLp, fJrelp, fJupp, cajsr, cansr, PhiCaL_ss, ICaL_ss, ICaL_i, I_ClCa, I_Clbk, ICaL_tot, Bcai, IKCa, Ta,d,ff,fs,fcaf,fcas,jca,nca];
end
end