#

# gotran file generated by cellml2gotran from ToRORd_dynCl_mid.cellml

parameters("Environment",
           celltype = 2)

parameters("Extracellular",
           nao = ScalarParam(140.0, unit="mM"),
           cao = ScalarParam(1.8, unit="mM"),
           ko = ScalarParam(5.0, unit="mM"),
           clo = ScalarParam(150.0, unit="mM"))

parameters("Physical constants",
           R = ScalarParam(8314.0, unit="J*kmole**-1*K**-1"),
           T = ScalarParam(310.0, unit="K"),
           F = ScalarParam(96485.0, unit="C*mole**-1"),
           zna = 1,
           zca = 2,
           zk = 1,
           zcl = -1)

parameters("Cell geometry",
           L = ScalarParam(0.01, unit="cm"),
           rad = ScalarParam(0.0011, unit="cm"))

states("CaMK",
       CaMKt = ScalarParam(0.0201882, unit="mM"))

parameters("CaMK",
           KmCaMK = ScalarParam(0.15, unit="mM"),
           aCaMK = ScalarParam(0.05, unit="mM**-1*ms**-1"),
           bCaMK = ScalarParam(0.00068, unit="ms**-1"),
           CaMKo = 0.05,
           KmCaM = ScalarParam(0.0015, unit="mM"))

parameters("Reversal potentials",
           PKNa = 0.01833)

parameters("IpCa",
           GpCa = ScalarParam(0.0005, unit="mS*uF**-1"),
           KmCap = ScalarParam(0.0005, unit="mM"))

parameters("Diff",
           tauNa = ScalarParam(2.0, unit="ms"),
           tauK = ScalarParam(2.0, unit="ms"),
           tauCa = ScalarParam(0.2, unit="ms"),
           tauCl = ScalarParam(2.0, unit="ms"))

parameters("INa", "Ca",
           INaCa_fractionSS = 0.35,
           kna1 = ScalarParam(15.0, unit="ms**-1"),
           kna2 = ScalarParam(5.0, unit="ms**-1"),
           kna3 = ScalarParam(88.12, unit="ms**-1"),
           kasymm = 12.5,
           wna = 6e4,
           wca = 6e4,
           wnaca = 5e3,
           kcaon = ScalarParam(1500000.0, unit="ms**-1"),
           kcaoff = ScalarParam(5000.0, unit="ms**-1"),
           qna = 0.5224,
           qca = 0.167,
           KmCaAct = ScalarParam(0.00015, unit="mM"),
           Gncx_b = ScalarParam(0.0034, unit="mS*uF**-1"))

parameters("INa", "K",
           k1p = ScalarParam(949.5, unit="ms**-1"),
           k1m = ScalarParam(182.4, unit="ms**-1"),
           k2p = ScalarParam(687.2, unit="ms**-1"),
           k2m = ScalarParam(39.4, unit="ms**-1"),
           k3p = ScalarParam(1899.0, unit="ms**-1"),
           k3m = ScalarParam(79300.0, unit="ms**-1"),
           k4p = ScalarParam(639.0, unit="ms**-1"),
           k4m = ScalarParam(40.0, unit="ms**-1"),
           Knai0 = ScalarParam(9.073, unit="mM"),
           Knao0 = ScalarParam(27.78, unit="mM"),
           delta = ScalarParam(-0.155, unit="mV"),
           Kki = ScalarParam(0.5, unit="ms**-1"),
           Kko = ScalarParam(0.3582, unit="ms**-1"),
           MgADP = ScalarParam(0.05, unit="mM"),
           MgATP = ScalarParam(9.8, unit="mM"),
           Kmgatp = ScalarParam(1.698e-07, unit="mM"),
           H = ScalarParam(1e-07, unit="mM"),
           eP = 4.2,
           Khp = ScalarParam(1.698e-07, unit="mM"),
           Knap = ScalarParam(224.0, unit="mM"),
           Kxkur = ScalarParam(292.0, unit="mM"),
           Pnak_b = ScalarParam(15.4509, unit="mS*uF**-1"))

parameters("INa", "b",
           PNab = ScalarParam(1.9239e-09, unit="mS*uF**-1"))

states("ICaL",
       d = -8.334604e-30,
       ff = 1.000000e+00,
       fs = 9.183587e-01,
       fcaf = 1.000000e+00,
       fcas = 9.997540e-01,
       jca = 9.999743e-01,
       ffp = 1.000000e+00,
       fcafp = 1.000000e+00,
       nca_ss = 5.336520e-04,
       nca_i = 1.257861e-03)

parameters("ICaL",
           Kmn = ScalarParam(0.002, unit="mM"),
           k2n = ScalarParam(500.0, unit="ms**-1"),
           PCa_b = 8.3757e-05,
           Aff = 0.6,
           tjca = ScalarParam(72.5, unit="ms"),
           vShift = ScalarParam(0.0, unit="mV"),
           offset = ScalarParam(0.0, unit="ms"),
           dielConstant = ScalarParam(74.0, unit="K**-1"),
           ICaL_fractionSS = 0.8)

parameters("SERCA",
           Jup_b = 1.0)

parameters("I_katp",
           gkatp = ScalarParam(4.3195, unit="mS*uF**-1"),
           fkatp = 0.0,
           K_o_n = ScalarParam(5.0, unit="mM"),
           A_atp = ScalarParam(2.0, unit="mM"),
           K_atp = ScalarParam(0.25, unit="mM"))

states("INa",
       m = 4.619565e-04,
       h = 8.739077e-01,
       j = 8.737841e-01,
       hp = 7.478972e-01,
       jp = 8.735375e-01)

parameters("INa",
           GNa = ScalarParam(11.7802, unit="mS*uF**-1"))

states("INa", "L",
       mL = 9.987709e-05,
       hL = 5.986118e-01,
       hLp = 3.339899e-01)

parameters("INa", "L",
           thL = ScalarParam(200.0, unit="ms"),
           GNaL_b = ScalarParam(0.0279, unit="mS*uF**-1"))

states("Ito",
       a = 7.994042e-04,
       iF = 9.997514e-01,
       iS = 5.702538e-01,
       ap = 4.072777e-04,
       iFp = 9.997514e-01,
       iSp = 6.351927e-01)

parameters("Ito",
           Gto_b = ScalarParam(0.16, unit="mS*uF**-1"),
           EKshift = ScalarParam(0.0, unit="mV"))

states("IKr",
       C3 = 5.910047e-04,
       C2 = 7.086460e-04,
       C1 = 9.983451e-01,
       O = 3.445733e-04,
       I = 1.064230e-05)

parameters("IKr",
           GKr_b = ScalarParam(0.0321, unit="mS*uF**-1"),
           alpha_1 = ScalarParam(0.154375, unit="ms**-1"),
           beta_1 = ScalarParam(0.1911, unit="ms**-1"))

states("IKs",
       xs1 = 2.642293e-01,
       xs2 = 1.327348e-04)

parameters("IKs",
           GKs_b = ScalarParam(0.0011, unit="mS*uF**-1"))

parameters("IK1",
           GK1_b = ScalarParam(0.6992, unit="mS*uF**-1"))

parameters("IKb",
           GKb_b = ScalarParam(0.0189, unit="mS*uF**-1"))

parameters("ICl",
           GClCa = ScalarParam(0.2843, unit="mS*uF**-1"),
           GClb = ScalarParam(0.00198, unit="mS*uF**-1"),
           KdClCa = ScalarParam(0.1, unit="mM"),
           Fjunc = 1)

parameters("ICab",
           PCab = ScalarParam(5.9194e-08, unit="mS*uF**-1"))

states("Ryr",
       Jrel_np = ScalarParam(-1.300486e-21, unit="mM*ms**-1"),
       Jrel_p = ScalarParam(-7.610714e-20, unit="mM*ms**-1"))

parameters("Ryr",
           bt = ScalarParam(4.75, unit="ms"),
           cajsr_half = ScalarParam(1.7, unit="mM"),
           Jrel_b = 1.5378)

states("Membrane",
       v = ScalarParam(-91.33918, unit="mV"))

parameters("Membrane",
           i_Stim_Start = ScalarParam(0.0, unit="ms"),
           i_Stim_End = ScalarParam(1e+17, unit="ms"),
           i_Stim_Amplitude = ScalarParam(-53.0, unit="uA*uF**-1"),
           i_Stim_Period = ScalarParam(1000.0, unit="ms"),
           i_Stim_PulseDuration = ScalarParam(1.0, unit="ms"))

parameters("Isac (Pueyo)--> ns + k",
Gsac_ns=0.006,
Esac_ns=-10,
lambda_max=1.1,
Gsac_k=0.2882*800/210)

expressions("Isac (Pueyo)--> ns + k")
Isac_P_ns=0
Isac_P_k=0


states("mechanics",
cai= ScalarParam(8.297576e-05, unit="mM"),
XS=0,
XW=0,
CaTrpn=0,
TmB=1,
Zetas=0,
Zetaw=0,
Cd=0)


states("Intracellular ions",
       nai = ScalarParam(15.94867, unit="mM"),
       nass = ScalarParam(15.94922, unit="mM"),
       ki = ScalarParam(156.7131, unit="mM"),
       kss = ScalarParam(156.713, unit="mM"),
       cli = ScalarParam(48.91277, unit="mM"),
       clss = ScalarParam(48.91274, unit="mM"),
       cass = ScalarParam(6.642187e-05, unit="mM"),
       cansr = ScalarParam(2.012225, unit="mM"),
       cajsr = ScalarParam(2.016415, unit="mM"))

parameters("Intracellular ions",
           cmdnmax_b = ScalarParam(0.05, unit="mM"),
           kmcmdn = ScalarParam(0.00238, unit="mM"),
           trpnmax = ScalarParam(0.07, unit="mM"),
           kmtrpn = ScalarParam(0.0005, unit="mM"),
           BSRmax = ScalarParam(0.047, unit="mM"),
           KmBSR = ScalarParam(0.00087, unit="mM"),
           BSLmax = ScalarParam(1.124, unit="mM"),
           KmBSL = ScalarParam(0.0087, unit="mM"),
           csqnmax = ScalarParam(10.0, unit="mM"),
           kmcsqn = ScalarParam(0.8, unit="mM"))

expressions("ToRORd dyn chloride")
vffrt = v*F*F/(R*T) # C*mole**-1
vfrt = v*F/(R*T)

expressions("Cell geometry")
vcell = 1000*3.14*rad*rad*L # ul
Ageo = 2*3.14*rad*rad + 2*3.14*rad*L # cm**2
Acap = 2*Ageo # cm**2
vmyo = 0.68*vcell # ul
vnsr = 0.0552*vcell # ul
vjsr = 0.0048*vcell # ul
vss = 0.02*vcell # ul

expressions("CaMK")
CaMKb = CaMKo*(1 - CaMKt)/(1 + KmCaM/cass) # mM
CaMKa = CaMKb + CaMKt # mM
dCaMKt_dt = aCaMK*CaMKb*(CaMKb + CaMKt) - bCaMK*CaMKt

expressions("Reversal potentials")
ENa = R*T/(zna*F)*log(nao/nai) # mV
EK = R*T/(zk*F)*log(ko/ki) # mV
EKs = R*T/(zk*F)*log((ko + PKNa*nao)/(ki + PKNa*nai)) # mV
ECl = R*T/(zcl*F)*log(clo/cli) # mV
EClss = R*T/(zcl*F)*log(clo/clss) # mV

expressions("IpCa")
IpCa = GpCa*cai/(KmCap + cai) # uA*uF**-1

expressions("Diff")
JdiffNa = (nass - nai)/tauNa # mM*ms**-1
JdiffK = (kss - ki)/tauK # mM*ms**-1
Jdiff = (cass - cai)/tauCa # mM*ms**-1
JdiffCl = (clss - cli)/tauNa # mM*ms**-1

expressions("Trans flux")
Jtr = (cansr - cajsr)/60 # mM*ms**-1

expressions("INa", "Ca")
hca = exp(qca*vfrt)
hna = exp(qna*vfrt)
h1_i = 1 + nai/kna3*(1 + hna)
h2_i = nai*hna/(kna3*h1_i)
h3_i = 1/h1_i
h4_i = 1 + nai/kna1*(1 + nai/kna2)
h5_i = nai*nai/(h4_i*kna1*kna2)
h6_i = 1/h4_i
h7_i = 1 + nao/kna3*(1 + 1/hna)
h8_i = nao/(kna3*hna*h7_i)
h9_i = 1/h7_i
h10_i = kasymm + 1 + nao/kna1*(1 + nao/kna2)
h11_i = nao*nao/(h10_i*kna1*kna2)
h12_i = 1/h10_i
k1_i = h12_i*cao*kcaon
k2_i = kcaoff
k3p_i = h9_i*wca
k3pp_i = h8_i*wnaca
k3_i = k3p_i + k3pp_i
k4p_i = h3_i*wca/hca
k4pp_i = h2_i*wnaca
k4_i = k4p_i + k4pp_i
k5_i = kcaoff
k6_i = h6_i*cai*kcaon
k7_i = h5_i*h2_i*wna
k8_i = h8_i*h11_i*wna
x1_i = k2_i*k4_i*(k7_i + k6_i) + k5_i*k7_i*(k2_i + k3_i)
x2_i = k1_i*k7_i*(k4_i + k5_i) + k4_i*k6_i*(k1_i + k8_i)
x3_i = k1_i*k3_i*(k7_i + k6_i) + k8_i*k6_i*(k2_i + k3_i)
x4_i = k2_i*k8_i*(k4_i + k5_i) + k3_i*k5_i*(k1_i + k8_i)
E1_i = x1_i/(x1_i + x2_i + x3_i + x4_i)
E2_i = x2_i/(x1_i + x2_i + x3_i + x4_i)
E3_i = x3_i/(x1_i + x2_i + x3_i + x4_i)
E4_i = x4_i/(x1_i + x2_i + x3_i + x4_i)
allo_i = 1/(1 + (KmCaAct/cai)**2)
JncxNa_i = 3*(E4_i*k7_i - E1_i*k8_i) + E3_i*k4pp_i - E2_i*k3pp_i # mM*ms**-1
JncxCa_i = E2_i*k2_i - E1_i*k1_i # mM*ms**-1
Gncx = Conditional(Eq(celltype, 1), Gncx_b*1.1, Conditional(Eq(celltype, 2), Gncx_b*1.4, Gncx_b)) # mS*uF**-1
INaCa_i = (1 - INaCa_fractionSS)*Gncx*allo_i*(zna*JncxNa_i + zca*JncxCa_i) # uA*uF**-1
h1_ss = 1 + nass/kna3*(1 + hna)
h2_ss = nass*hna/(kna3*h1_ss)
h3_ss = 1/h1_ss
h4_ss = 1 + nass/kna1*(1 + nass/kna2)
h5_ss = nass*nass/(h4_ss*kna1*kna2)
h6_ss = 1/h4_ss
h7_ss = 1 + nao/kna3*(1 + 1/hna)
h8_ss = nao/(kna3*hna*h7_ss)
h9_ss = 1/h7_ss
h10_ss = kasymm + 1 + nao/kna1*(1 + nao/kna2)
h11_ss = nao*nao/(h10_ss*kna1*kna2)
h12_ss = 1/h10_ss
k1_ss = h12_ss*cao*kcaon
k2_ss = kcaoff
k3p_ss = h9_ss*wca
k3pp_ss = h8_ss*wnaca
k3_ss = k3p_ss + k3pp_ss
k4p_ss = h3_ss*wca/hca
k4pp_ss = h2_ss*wnaca
k4_ss = k4p_ss + k4pp_ss
k5_ss = kcaoff
k6_ss = h6_ss*cass*kcaon
k7_ss = h5_ss*h2_ss*wna
k8_ss = h8_ss*h11_ss*wna
x1_ss = k2_ss*k4_ss*(k7_ss + k6_ss) + k5_ss*k7_ss*(k2_ss + k3_ss)
x2_ss = k1_ss*k7_ss*(k4_ss + k5_ss) + k4_ss*k6_ss*(k1_ss + k8_ss)
x3_ss = k1_ss*k3_ss*(k7_ss + k6_ss) + k8_ss*k6_ss*(k2_ss + k3_ss)
x4_ss = k2_ss*k8_ss*(k4_ss + k5_ss) + k3_ss*k5_ss*(k1_ss + k8_ss)
E1_ss = x1_ss/(x1_ss + x2_ss + x3_ss + x4_ss)
E2_ss = x2_ss/(x1_ss + x2_ss + x3_ss + x4_ss)
E3_ss = x3_ss/(x1_ss + x2_ss + x3_ss + x4_ss)
E4_ss = x4_ss/(x1_ss + x2_ss + x3_ss + x4_ss)
allo_ss = 1/(1 + (KmCaAct/cass)**2)
JncxNa_ss = 3*(E4_ss*k7_ss - E1_ss*k8_ss) + E3_ss*k4pp_ss - E2_ss*k3pp_ss # mM*ms**-1
JncxCa_ss = E2_ss*k2_ss - E1_ss*k1_ss # mM*ms**-1
INaCa_ss = INaCa_fractionSS*Gncx*allo_ss*(zna*JncxNa_ss + zca*JncxCa_ss) # uA*uF**-1

expressions("INa", "K")
Knai = Knai0*exp(delta*vfrt/3) # mM
Knao = Knao0*exp((1 - delta)*vfrt/3) # mM
P = eP/(1 + H/Khp + nai/Knap + ki/Kxkur)
a1 = k1p*(nai/Knai)**3/((1 + nai/Knai)**3 + (1 + ki/Kki)**2 - 1)
b1 = k1m*MgADP
a2 = k2p
b2 = k2m*(nao/Knao)**3/((1 + nao/Knao)**3 + (1 + ko/Kko)**2 - 1)
a3 = k3p*(ko/Kko)**2/((1 + nao/Knao)**3 + (1 + ko/Kko)**2 - 1)
b3 = k3m*P*H/(1 + MgATP/Kmgatp)
a4 = k4p*MgATP/Kmgatp/(1 + MgATP/Kmgatp)
b4 = k4m*(ki/Kki)**2/((1 + nai/Knai)**3 + (1 + ki/Kki)**2 - 1)
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*(E1*a3 - E2*b3) # mM*ms**-1
JnakK = 2*(E4*b1 - E3*a1) # mM*ms**-1
Pnak = Conditional(Eq(celltype, 1), Pnak_b*0.9, Conditional(Eq(celltype, 2), Pnak_b*0.7, Pnak_b)) # mS*uF**-1
INaK = Pnak*(zna*JnakNa + zk*JnakK) # uA*uF**-1

expressions("INa", "b")
INab = PNab*vffrt*(nai*exp(vfrt) - nao)/(exp(vfrt) - 1) # uA*uF**-1

expressions("ICaL")
dss = Conditional(Ge(v, 31.4978), 1, 1.0763*exp(-1.0070*exp(-0.0829*v)))
td = offset + 0.6 + 1/(exp(-0.05*(v + vShift + 6)) + exp(0.09*(v + vShift + 14))) # ms
dd_dt = (dss - d)/td
fss = 1/(1 + exp((v + 19.58)/3.696))
tff = 7 + 1/(0.0045*exp(-(v + 20)/10) + 0.0045*exp((v + 20)/10)) # ms
tfs = 1000 + 1/(0.000035*exp(-(v + 5)/4) + 0.000035*exp((v + 5)/6)) # ms
Afs = 1 - Aff
dff_dt = (fss - ff)/tff
dfs_dt = (fss - fs)/tfs
f = Aff*ff + Afs*fs
fcass = fss
tfcaf = 7 + 1/(0.04*exp(-(v - 4)/7) + 0.04*exp((v - 4)/7)) # ms
tfcas = 100 + 1/(0.00012*exp(-v/3) + 0.00012*exp(v/7)) # ms
Afcaf = 0.3 + 0.6/(1 + exp((v - 10)/10))
Afcas = 1 - Afcaf
dfcaf_dt = (fcass - fcaf)/tfcaf
dfcas_dt = (fcass - fcas)/tfcas
fca = Afcaf*fcaf + Afcas*fcas
jcass = 1.0/(1.0 + exp((v + 18.08)/2.7916))
djca_dt = (jcass - jca)/tjca
tffp = 2.5*tff # ms
dffp_dt = (fss - ffp)/tffp
fp = Aff*ffp + Afs*fs
tfcafp = 2.5*tfcaf # ms
dfcafp_dt = (fcass - fcafp)/tfcafp
fcap = Afcaf*fcafp + Afcas*fcas
km2n = jca*1 # ms**-1
anca_ss = 1/(k2n/km2n + (1 + Kmn/cass)**4)
dnca_ss_dt = anca_ss*k2n - nca_ss*km2n
Io = 0.5*(nao + ko + clo + 4*cao)/1000
Iss = 0.5*(nass + kss + clss + 4*cass)/1000
constA = 1.82e6*(dielConstant*T)**-1.5
gamma_cass = exp(-constA*4*(sqrt(Iss)/(1 + sqrt(Iss)) - 0.3*Iss))
gamma_cao = exp(-constA*4*(sqrt(Io)/(1 + sqrt(Io)) - 0.3*Io))
gamma_nass = exp(-constA*1*(sqrt(Iss)/(1 + sqrt(Iss)) - 0.3*Iss))
gamma_nao = exp(-constA*1*(sqrt(Io)/(1 + sqrt(Io)) - 0.3*Io))
gamma_kss = exp(-constA*1*(sqrt(Iss)/(1 + sqrt(Iss)) - 0.3*Iss))
gamma_ko = exp(-constA*1*(sqrt(Io)/(1 + sqrt(Io)) - 0.3*Io))
PhiCaL_ss = 4*vffrt*(gamma_cass*cass*exp(2*vfrt) - gamma_cao*cao)/(exp(2*vfrt) - 1)
PhiCaNa_ss = 1*vffrt*(gamma_nass*nass*exp(1*vfrt) - gamma_nao*nao)/(exp(1*vfrt) - 1)
PhiCaK_ss = 1*vffrt*(gamma_kss*kss*exp(1*vfrt) - gamma_ko*ko)/(exp(1*vfrt) - 1)
PCa = Conditional(Eq(celltype, 1), PCa_b*1.2, Conditional(Eq(celltype, 2), PCa_b*2, PCa_b))
PCap = 1.1*PCa
PCaNa = 0.00125*PCa
PCaK = 3.574e-4*PCa
PCaNap = 0.00125*PCap
PCaKp = 3.574e-4*PCap
fICaLp = 1/(1 + KmCaMK/CaMKa)
ICaL_ss = ICaL_fractionSS*((1 - fICaLp)*PCa*PhiCaL_ss*d*(f*(1 - nca_ss) + jca*fca*nca_ss) + fICaLp*PCap*PhiCaL_ss*d*(fp*(1 - nca_ss) + jca*fcap*nca_ss)) # uA*uF**-1
ICaNa_ss = ICaL_fractionSS*((1 - fICaLp)*PCaNa*PhiCaNa_ss*d*(f*(1 - nca_ss) + jca*fca*nca_ss) + fICaLp*PCaNap*PhiCaNa_ss*d*(fp*(1 - nca_ss) + jca*fcap*nca_ss)) # uA*uF**-1
ICaK_ss = ICaL_fractionSS*((1 - fICaLp)*PCaK*PhiCaK_ss*d*(f*(1 - nca_ss) + jca*fca*nca_ss) + fICaLp*PCaKp*PhiCaK_ss*d*(fp*(1 - nca_ss) + jca*fcap*nca_ss)) # uA*uF**-1
anca_i = 1/(k2n/km2n + (1 + Kmn/cai)**4)
dnca_i_dt = anca_i*k2n - nca_i*km2n
Ii = 0.5*(nai + ki + cli + 4*cai)/1000
gamma_cai = exp(-constA*4*(sqrt(Ii)/(1 + sqrt(Ii)) - 0.3*Ii))
gamma_nai = exp(-constA*1*(sqrt(Ii)/(1 + sqrt(Ii)) - 0.3*Ii))
gamma_ki = exp(-constA*1*(sqrt(Ii)/(1 + sqrt(Ii)) - 0.3*Ii))
PhiCaL_i = 4*vffrt*(gamma_cai*cai*exp(2*vfrt) - gamma_cao*cao)/(exp(2*vfrt) - 1)
PhiCaNa_i = 1*vffrt*(gamma_nai*nai*exp(1*vfrt) - gamma_nao*nao)/(exp(1*vfrt) - 1)
PhiCaK_i = 1*vffrt*(gamma_ki*ki*exp(1*vfrt) - gamma_ko*ko)/(exp(1*vfrt) - 1)
ICaL_i = (1 - ICaL_fractionSS)*((1 - fICaLp)*PCa*PhiCaL_i*d*(f*(1 - nca_i) + jca*fca*nca_i) + fICaLp*PCap*PhiCaL_i*d*(fp*(1 - nca_i) + jca*fcap*nca_i)) # uA*uF**-1
ICaNa_i = (1 - ICaL_fractionSS)*((1 - fICaLp)*PCaNa*PhiCaNa_i*d*(f*(1 - nca_i) + jca*fca*nca_i) + fICaLp*PCaNap*PhiCaNa_i*d*(fp*(1 - nca_i) + jca*fcap*nca_i)) # uA*uF**-1
ICaK_i = (1 - ICaL_fractionSS)*((1 - fICaLp)*PCaK*PhiCaK_i*d*(f*(1 - nca_i) + jca*fca*nca_i) + fICaLp*PCaKp*PhiCaK_i*d*(fp*(1 - nca_i) + jca*fcap*nca_i)) # uA*uF**-1
ICaL = ICaL_ss + ICaL_i # uA*uF**-1
ICaNa = ICaNa_ss + ICaNa_i # uA*uF**-1
ICaK = ICaK_ss + ICaK_i # uA*uF**-1

expressions("SERCA")
upScale = Conditional(Eq(celltype, 1), 1.3, 1)
Jupnp = upScale*0.005425*cai/(cai + 0.00092) # mM*ms**-1
Jupp = upScale*2.75*0.005425*cai/(cai + 0.00092 - 0.00017) # mM*ms**-1
fJupp = 1/(1 + KmCaMK/CaMKa)
Jleak = 0.0048825*cansr/15 # mM*ms**-1
Jup = Jup_b*((1 - fJupp)*Jupnp + fJupp*Jupp - Jleak) # mM*ms**-1

expressions("I_katp")
akik = (ko/K_o_n)**0.24
bkik = 1/(1 + (A_atp/K_atp)**2)
I_katp = fkatp*gkatp*akik*bkik*(v - EK) # uA*uF**-1

expressions("INa")
mss = 1/(1 + exp(-(v + 56.86)/9.03))**2
tm = 0.1292*exp(-((v + 45.79)/15.54)**2) + 0.06487*exp(-((v - 4.823)/51.12)**2) # ms
dm_dt = (mss - m)/tm
hss = 1/(1 + exp((v + 71.55)/7.43))**2
ah = Conditional(Ge(v, -40), 0, 0.057*exp(-(v + 80)/6.8))
bh = Conditional(Ge(v, -40), 0.77/(0.13*(1 + exp(-(v + 10.66)/11.1))), 2.7*exp(0.079*v) + 3.1e5*exp(0.3485*v))
th = 1/(ah + bh) # ms
dh_dt = (hss - h)/th
aj = Conditional(Ge(v, -40), 0, (-2.5428e4*exp(0.2444*v) - 6.948e-6*exp(-0.04391*v))*(v + 37.78)/(1 + exp(0.311*(v + 79.23))))
bj = Conditional(Ge(v, -40), 0.6*exp(0.057*v)/(1 + exp(-0.1*(v + 32))), 0.02424*exp(-0.01052*v)/(1 + exp(-0.1378*(v + 40.14))))
jss = hss
tj = 1/(aj + bj) # ms
dj_dt = (jss - j)/tj
hssp = 1/(1 + exp((v + 77.55)/7.43))**2
dhp_dt = (hssp - hp)/th
tjp = 1.46*tj # ms
djp_dt = (jss - jp)/tjp
fINap = 1/(1 + KmCaMK/CaMKa)
INa = GNa*(v - ENa)*m**3*((1 - fINap)*h*j + fINap*hp*jp) # uA*uF**-1

expressions("INa", "L")
mLss = 1/(1 + exp(-(v + 42.85)/5.264))
tmL = 0.1292*exp(-((v + 45.79)/15.54)**2) + 0.06487*exp(-((v - 4.823)/51.12)**2) # ms
dmL_dt = (mLss - mL)/tmL
hLss = 1/(1 + exp((v + 87.61)/7.488))
dhL_dt = (hLss - hL)/thL
hLssp = 1/(1 + exp((v + 93.81)/7.488))
thLp = 3*thL # ms
dhLp_dt = (hLssp - hLp)/thLp
GNaL = Conditional(Eq(celltype, 1), GNaL_b*0.6, GNaL_b) # mS*uF**-1
fINaLp = 1/(1 + KmCaMK/CaMKa)
INaL = GNaL*(v - ENa)*mL*((1 - fINaLp)*hL + fINaLp*hLp) # uA*uF**-1

expressions("Ito")
ass = 1/(1 + exp(-(v + EKshift - 14.34)/14.82))
ta = 1.0515/(1/(1.2089*(1 + exp(-(v + EKshift - 18.4099)/29.3814))) + 3.5/(1 + exp((v + EKshift + 100)/29.3814))) # ms
da_dt = (ass - a)/ta
iss = 1/(1 + exp((v + EKshift + 43.94)/5.711))
delta_epi = Conditional(Eq(celltype, 1), 1 - 0.95/(1 + exp((v + EKshift + 70)/5)), 1)
tiF_b = 4.562 + 1/(0.3933*exp(-(v + EKshift + 100)/100) + 0.08004*exp((v + EKshift + 50)/16.59)) # ms
tiS_b = 23.62 + 1/(0.001416*exp(-(v + EKshift + 96.52)/59.05) + 1.78e-8*exp((v + EKshift + 114.1)/8.079)) # ms
tiF = tiF_b*delta_epi # ms
tiS = tiS_b*delta_epi # ms
AiF = 1/(1 + exp((v + EKshift - 213.6)/151.2))
AiS = 1 - AiF
diF_dt = (iss - iF)/tiF
diS_dt = (iss - iS)/tiS
i = AiF*iF + AiS*iS
assp = 1/(1 + exp(-(v + EKshift - 24.34)/14.82))
dap_dt = (assp - ap)/ta
dti_develop = 1.354 + 1e-4/(exp((v + EKshift - 167.4)/15.89) + exp(-(v + EKshift - 12.23)/0.2154))
dti_recover = 1 - 0.5/(1 + exp((v + EKshift + 70.0)/20.0))
tiFp = dti_develop*dti_recover*tiF # ms
tiSp = dti_develop*dti_recover*tiS # ms
diFp_dt = (iss - iFp)/tiFp
diSp_dt = (iss - iSp)/tiSp
ip = AiF*iFp + AiS*iSp
Gto = Conditional(Eq(celltype, 1), Gto_b*2, Conditional(Eq(celltype, 2), Gto_b*2, Gto_b)) # mS*uF**-1
fItop = 1/(1 + KmCaMK/CaMKa)
Ito = Gto*(v - EK)*((1 - fItop)*a*i + fItop*ap*ip) # uA*uF**-1

expressions("IKr")
alpha = 0.1161*exp(0.2990*vfrt) # ms**-1
beta = 0.2442*exp(-1.604*vfrt) # ms**-1
alpha_2 = 0.0578*exp(0.9710*vfrt) # ms**-1
beta_2 = 0.349e-3*exp(-1.062*vfrt) # ms**-1
alpha_i = 0.2533*exp(0.5953*vfrt) # ms**-1
beta_i = 0.06525*exp(-0.8209*vfrt) # ms**-1
alpha_C2ToI = 0.52e-4*exp(1.525*vfrt) # ms**-1
beta_ItoC2 = beta_2*beta_i*alpha_C2ToI/(alpha_2*alpha_i) # ms**-1
dC3_dt = beta*C2 - alpha*C3
dC2_dt = alpha*C3 + beta_1*C1 - (beta + alpha_1)*C2
dC1_dt = alpha_1*C2 + beta_2*O + beta_ItoC2*I - (beta_1 + alpha_2 + alpha_C2ToI)*C1
dO_dt = alpha_2*C1 + beta_i*I - (beta_2 + alpha_i)*O
dI_dt = alpha_C2ToI*C1 + alpha_i*O - (beta_ItoC2 + beta_i)*I
GKr = Conditional(Eq(celltype, 1), GKr_b*1.3, Conditional(Eq(celltype, 2), GKr_b*0.8, GKr_b)) # mS*uF**-1
IKr = GKr*sqrt(ko/5)*O*(v - EK) # uA*uF**-1

expressions("IKs")
xs1ss = 1/(1 + exp(-(v + 11.6)/8.932))
txs1 = 817.3 + 1/(2.326e-4*exp((v + 48.28)/17.8) + 0.001292*exp(-(v + 210)/230)) # ms
dxs1_dt = (xs1ss - xs1)/txs1
xs2ss = xs1ss
txs2 = 1/(0.01*exp((v - 50)/20) + 0.0193*exp(-(v + 66.54)/31)) # ms
dxs2_dt = (xs2ss - xs2)/txs2
KsCa = 1 + 0.6/(1 + (3.8e-5/cai)**1.4)
GKs = Conditional(Eq(celltype, 1), GKs_b*1.4, GKs_b) # mS*uF**-1
IKs = GKs*KsCa*xs1*xs2*(v - EKs) # uA*uF**-1

expressions("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 = Conditional(Eq(celltype, 1), GK1_b*1.2, Conditional(Eq(celltype, 2), GK1_b*1.3, GK1_b)) # mS*uF**-1
IK1 = GK1*sqrt(ko/5)*K1ss*(v - EK) # uA*uF**-1

expressions("IKb")
xkb = 1/(1 + exp(-(v - 10.8968)/23.9871))
GKb = Conditional(Eq(celltype, 1), GKb_b*0.6, GKb_b) # mS*uF**-1
IKb = GKb*xkb*(v - EK) # uA*uF**-1

expressions("ICl")
IClCa_junc = Fjunc*GClCa/(1 + KdClCa/cass)*(v - EClss) # uA*uF**-1
IClCa_sl = (1 - Fjunc)*GClCa/(1 + KdClCa/cai)*(v - ECl) # uA*uF**-1
IClCa = IClCa_junc + IClCa_sl # uA*uF**-1
IClb = GClb*(v - ECl) # uA*uF**-1

expressions("ICab")
ICab = PCab*4*vffrt*(gamma_cai*cai*exp(2*vfrt) - gamma_cao*cao)/(exp(2*vfrt) - 1) # uA*uF**-1

expressions("Ryr")
a_rel = 0.5*bt/1 # mM*ms**-1
Jrel_inf_b = -a_rel*ICaL_ss/1/(1 + (cajsr_half/cajsr)**8) # mM*ms**-1
Jrel_inf = Conditional(Eq(celltype, 2), Jrel_inf_b*1.7, Jrel_inf_b) # mM*ms**-1
tau_rel_b = bt/(1 + 0.0123/cajsr) # ms
tau_rel = Conditional(Lt(tau_rel_b, 0.001), 0.001, tau_rel_b) # ms
dJrel_np_dt = (Jrel_inf - Jrel_np)/tau_rel
btp = 1.25*bt # ms
a_relp = 0.5*btp/1 # mM*ms**-1
Jrel_infp_b = -a_relp*ICaL_ss/1/(1 + (cajsr_half/cajsr)**8) # mM*ms**-1
Jrel_infp = Conditional(Eq(celltype, 2), Jrel_infp_b*1.7, Jrel_infp_b) # mM*ms**-1
tau_relp_b = btp/(1 + 0.0123/cajsr) # ms
tau_relp = Conditional(Lt(tau_relp_b, 0.001), 0.001, tau_relp_b) # ms
dJrel_p_dt = (Jrel_infp - Jrel_p)/tau_relp
fJrelp = 1/(1 + KmCaMK/CaMKa)
Jrel = Jrel_b*((1 - fJrelp)*Jrel_np + fJrelp*Jrel_p) # mM*ms**-1

expressions("Membrane")
Istim = Conditional(And(Ge(time, i_Stim_Start), Le(time - i_Stim_Start - floor((time - i_Stim_Start)/i_Stim_Period)*i_Stim_Period, i_Stim_PulseDuration), ), i_Stim_Amplitude, 0) # uA*uF**-1
dv_dt = -(INa + INaL + Ito + ICaL + ICaNa + ICaK + IKr + IKs + IK1 + INaCa_i + INaCa_ss + INaK + INab + IKb + IpCa + ICab + IClCa + IClb + I_katp + Istim + Isac_P_ns + Isac_P_k)

expressions("Intracellular ions")
cmdnmax = Conditional(Eq(celltype, 1), cmdnmax_b*1.3, cmdnmax_b) # mM
dnai_dt = -(INa + INaL + 3*INaCa_i + ICaNa_i + 3*INaK + INab + Isac_P_ns/3)*Acap/(F*vmyo) + JdiffNa*vss/vmyo
dnass_dt = -(ICaNa_ss + 3*INaCa_ss)*Acap/(F*vss) - JdiffNa
dki_dt = -(Ito + IKr + IKs + IK1 + IKb + I_katp + Istim - 2*INaK + ICaK_i + Isac_P_k+Isac_P_ns/3)*Acap/(F*vmyo) + JdiffK*vss/vmyo
dkss_dt = -ICaK_ss*Acap/(F*vss) - JdiffK
dcli_dt = (IClb + IClCa_sl)*Acap/(F*vmyo) + JdiffCl*vss/vmyo
dclss_dt = IClCa_junc*Acap/(F*vss) - JdiffCl
Bcai = 1/(1 + cmdnmax*kmcmdn/(kmcmdn + cai)**2 + trpnmax*kmtrpn/(kmtrpn + cai)**2)

Bcass = 1/(1 + BSRmax*KmBSR/(KmBSR + cass)**2 + BSLmax*KmBSL/(KmBSL + cass)**2)
dcass_dt = Bcass*(-(ICaL_ss - 2*INaCa_ss)*Acap/(2*F*vss) + Jrel*vjsr/vss - Jdiff)
dcansr_dt = Jup - Jtr*vjsr/vnsr
Bcajsr = 1/(1 + csqnmax*kmcsqn/(kmcsqn + cajsr)**2)
dcajsr_dt = Bcajsr*(Jtr - Jrel)


parameters("mechanics",
emcoupling=1,
lmbda=1,
dLambda=0,
mode=1,
isacs=0,
calib=1,
ktrpn = 0.1,
ntrpn = 2,
Trpn50 = 0.35,
rw = 0.5,
rs = 0.25,
gammas = 0.0085,
gammaw = 0.615,
phi = 2.23,
Tot_A = 25,
Beta0 = 2.3,
Beta1 = -2.4,
cat50_ref = 0.805,
Tref = 120,
kuw = 0.182,
kws = 0.012,
ku=0.04,
ntm=2.4,
p_a = 2.1,
p_b = 9.1,
p_k = 7,
etal = 200,
etas = 20)


expressions("mechanics")
XS_max = Conditional(Gt(XS, 0), XS, 0)
XW_max = Conditional(Gt(XW, 0), XW, 0)
CaTrpn_max = Conditional(Gt(CaTrpn, 0), CaTrpn, 0)
kwu = kuw*(1/rw-1)-kws
ksu = kws*rw*(1/rs-1)
Aw = Tot_A*rs/((1-rs)*rw+rs)
As = Aw
cw = phi*kuw*((1-rs)*(1-rw))/((1-rs)*rw)
cs = phi*kws*((1-rs)*rw)/rs
lambda_min12 = Conditional(Lt(lmbda, 1.2), lmbda, 1.2)
lambda_min087 = Conditional(Lt(lambda_min12, 0.87), lambda_min12, 0.87)
h_lambda_prima = 1+Beta0*(lambda_min12+lambda_min087-1.87)
h_lambda = Conditional(Gt(h_lambda_prima, 0), h_lambda_prima, 0)
XU = (1-TmB)-XS-XW
gammawu = gammaw*abs(Zetaw)
gammasu  = gammas*Conditional(Gt(Gt(Zetas,0)*Zetas, Lt(Zetas, -1)*(-Zetas-1)), Gt(Zetas,0)*Zetas, Lt(Zetas, -1)*(-Zetas-1))
dXS_dt = kws*XW - ksu*XS - gammasu*XS
dXW_dt = kuw*XU - kwu*XW - kws*XW - gammawu*XW
cat50 = (cat50_ref+Beta1*(lambda_min12-1))
dCaTrpn_dt = ktrpn*(((cai*1000/cat50)**ntrpn)*(1-CaTrpn)-CaTrpn)
kb = ku*Trpn50**ntm/(1-rs-(1-rs)*rw)
dTmB_dt = kb*Conditional(Lt(CaTrpn**(-ntm/2), 100), CaTrpn**(-ntm/2), 100)*XU-ku*CaTrpn**(ntm/2)*TmB
dZetas_dt = As*dLambda - cs*Zetas
dZetaw_dt = Aw*dLambda - cw*Zetaw
Ta = h_lambda*(Tref/rs)*(XS*(Zetas+1) + XW*Zetaw)
C = lambda_min12 - 1
dCd = C - Cd
eta = Conditional(Lt(dCd, 0), etas, etal)
dCd_dt = p_k * (C - Cd) / eta
Fd = eta * dCd
F1 = (exp( p_b * C) - 1)
Tp = p_a * (F1 + Fd)
Ttot = Ta + Tp
J_TRPN = dCaTrpn_dt*trpnmax
dcai_dt=Bcai*(-(IpCa+ICab-2.0*INaCa_i+Isac_P_ns/3)*Acap/(2.0*F*vmyo)-Jup*vnsr/vmyo+Jdiff*vss/vmyo - J_TRPN)