ellinwood-2017.mmt

[[model]]
name: ellinwood-2017
version: 20240904
mmt_authors: Michael Clerx
display_name: Ellinwood et al., 2017
desc: """
    The 2017 model of the human atrial AP by Ellinwood et al. [1, 2].
    Based on the adaptation of the Grandi model [3] by Morotti et al. [4].

    This implementation is based on the original Matlab code provided on the
    Grandi lab website [5]. It was checked against the original code numerically
    by comparing the calculated derivatives.

    The stimulus was set to 0.5 ms and approximately twice the threshold value
    for depolarisation.

    References:

    [1] Ellinwood, Dobrev, Morotti, Grandi (2017) Revealing kinetics and
        state-dependent binding properties of IKur-targeting drugs that
        maximize atrial fibrillation selectivity. Chaos, 27, 093918
        https://doi.org/10.1063/1.5000226

    [2] Erratum to [1]
        https://doi.org/10.1063/1.5007051

    [3] *Grandi, E., *Pandit, S. V., *Voigt, N., Workman, A. J., Dobrev, D.,
        Jalife, J., & Bers, D. M. (2011). *Shared first authors. Human Atrial
        Action Potential and Ca2+ Model: Sinus Rhythm and Chronic Atrial
        Fibrillation. Circulation Research, 109(9), 1055-1066.
        https://doi.org/10.1161/CIRCRESAHA.111.253955

    [4] Morotti, S., McCulloch, A. D., Bers, D. M., Edwards, A. G.,
        & Grandi, E. (2016) Journal of Molecular and Cellular Cardiology
        https://doi.org/10.1016/j.yjmcc.2015.07.030

    [5] /~https://github.com/drgrandilab/Ellinwood-et-al-2017-Human-Atrial-Model-with-Updated-IKur
        Accessed on 2024-08-07

    """
# Initial values
membrane.V    = -7.41648369976224160e+01
ina.IC3       =  1.52061911704855252e-01
ina.IC2       =  1.21792676095848484e-02
ina.IF        =  4.54740797137192111e-04
ina.IM1       =  4.83257687571400363e-01
ina.C3        =  3.24316241134374128e-01
ina.C2        =  2.58882106924117616e-02
ina.C1        =  8.96316106602432694e-04
ina.O         =  6.82949304728142615e-06
ina.LC3       =  3.07947805371166006e-07
ina.LC2       =  2.45814664140850826e-08
ina.LC1       =  8.51017457460864658e-10
ina.LO        =  6.47009541610394701e-12
ical.d        =  1.91966101244576468e-05
ical.f        =  9.98579106751354639e-01
ical.fCa_jn   =  4.72106407222124685e-02
ical.fCa_sl   =  3.40031498654320141e-02
ito.x         =  1.29068265624421460e-03
ito.y         =  9.49128113097513548e-01
ikr.x         =  1.88822070115933428e-03
iks.x         =  7.11989787652290550e-03
ikur.C1       =  3.34704616007602318e-01
ikur.C2       =  4.08158462163515590e-01
ikur.C3       =  1.86652599711293971e-01
ikur.C4       =  3.79396227844088704e-02
ikur.O        =  2.89484795588510098e-03
ryr.r         =  8.01184984859831362e-01
ryr.o         =  2.35128842890773916e-06
ryr.i         =  5.83493121959741624e-07
nabss.NaB_jn  =  3.70554209310082383e+00
nabss.NaB_sl  =  8.08608652600110056e-01
cabct.TnCL    =  1.88259407185373505e-02
cabct.TnCHCa  =  1.28707878453568375e-01
cabct.TnCHMg  =  5.28487360732295911e-03
cabct.CaM     =  7.28983716096573853e-04
cabct.MyoCa   =  4.25383768990801776e-03
cabct.MyoMg   =  1.35230335919070938e-01
cabct.SRB     =  4.58091058990180676e-03
cabss.SLL_jn  =  1.42915067546659797e-02
cabss.SLL_sl  =  2.22278011276465619e-02
cabss.SLH_jn  =  1.06310545025060987e-01
cabss.SLH_sl  =  1.94815512517596906e-01
cabsr.Csqn    =  1.15780761298856683e+00
cabsr.Ca_sr   =  5.21827747395823160e-01
sodium.Na_jn  =  9.60544045288605730e+00
sodium.Na_sl  =  9.60598913412345468e+00
sodium.Na_i   =  9.60619407675032022e+00
calcium.Ca_jn =  3.43756620674409621e-04
calcium.Ca_sl =  2.42132268244438914e-04
calcium.Ca_i  =  2.19152503740116723e-04

#
# Simulation
#
[engine]
time = 0 [ms]
    in [ms]
    bind time
pace = 0
    bind pace

#
# Model modes
#
[mode]
AF = 0
    desc: Atrial fibrillation remodelling yes/no
ISO = 0
    desc: Apply isoproterenol yes/no
ACh = 0 [mM] in [mM]
    desc: Acetylcholine concentration

#
# Membrane potential
#
[membrane]
use stimulus.i_stim, potassium.I_K_tot
dot(V) = -(i_ion + i_stim + i_diff)
    in [mV]
    label membrane_potential
I_Na_tot = sodium.I_Na_tot_jn + sodium.I_Na_tot_sl
    in [A/F]
I_Ca_tot = calcium.I_Ca_tot_jn + calcium.I_Ca_tot_sl
    in [A/F]
I_Cl_tot = iclca.IClCa + iclb.IClB
    in [A/F]
i_ion = I_Na_tot + I_Cl_tot + I_Ca_tot + I_K_tot
    in [A/F]
    label cellular_current
i_diff = 0 [A/F]
    in [A/F]
    bind diffusion_current

#
# Stimulus current
#
[stimulus]
i_stim = engine.pace * amplitude
    in [A/F]
    label stimulus_current
amplitude = -19 [A/F] * 2
    in [A/F]

#
# Physical constants
#
[phys]
F = 96485 [C/mol]
    in [C/mol]
R = 8314 [mJ/mol/K]
    in [mJ/mol/K]
T = 310 [K]
    in [K]
RTF = R * T / F
    in [mV]
FRT = F / R / T
    in [1/mV]
Q = (T - 310 [K]) / 10 [K]

#
# Cell geometry
#
[cell]
C = 1.1e-10 [F]
    in [F]
    label membrane_capacitance
L = 100 [um]
    in [um]
R = 10.25 [um]
    in [um]
pi = 3.141592653589793
Vcell = pi * R^2 * L * 1e-15 [L/um^3]
    in [L]
Vmyo = 0.65 * Vcell
    in [L]
Vsr = 0.035 * Vcell
    in [L]
Vsl = 0.02 * Vcell
    in [L]
Vjn = 0.0539 * 0.01 * Vcell
    in [L]
fjn = 0.11
    desc: Fraction of currents in junctional compartment
fsl = 1 - fjn
    desc: Fraction of currents in subsarcolemmal compartment
fjn_CaL = 0.9
    desc: Fraction of ICaL in junctional compartment
fsl_CaL = 1 - fjn_CaL
    desc: Fraction of ICaL in subsarcolemmal compartment

#
# Fixed ion concentrations
#
[ion]
Na_o = 140 [mM]
    desc: Extracellular Na
    in [mM]
K_o = 5.4 [mM]
    desc: Extracellular K
    in [mM]
    label K_o
K_i = 120 [mM]
    desc: Intracellular K
    in [mM]
Ca_o = 1.8 [mM]
    desc: Extracellular Ca
    in [mM]
Cl_i = 15 [mM]
    desc: Intracellular Cl
    in [mM]
Cl_o = 150 [mM]
    desc: Extracellular Cl
    in [mM]
Mg_i = 1 [mM]
    desc: Intracellular Mg
    in [mM]

#
# Reversal potentials
#
[rev]
use phys.RTF
ENa_jn = RTF * log(ion.Na_o / sodium.Na_jn)
    in [mV]
ENa_sl = RTF * log(ion.Na_o / sodium.Na_sl)
    in [mV]
ECa_jn = 0.5 * RTF * log(ion.Ca_o / calcium.Ca_jn)
    in [mV]
ECa_sl = 0.5 * RTF * log(ion.Ca_o / calcium.Ca_sl)
    in [mV]
ECl = RTF * log(ion.Cl_i / ion.Cl_o)
    in [mV]
EK = RTF * log(ion.K_o / ion.K_i)
    in [mV]
EKs = RTF * log((ion.K_o + pNaK * ion.Na_o) / (ion.K_i + pNaK * sodium.Na_i))
    in [mV]
pNaK = 0.01833

#
# Fast sodium current (drug-free version)
#
# Added in Morotti 2016.
#
[ina]
use membrane.V
INa = INa_jn + INa_sl
    in [A/F]
INa_jn = cell.fjn * gNa * (O + LO) * (V - rev.ENa_jn)
    in [A/F]
INa_sl = cell.fsl * gNa * (O + LO) * (V - rev.ENa_sl)
    in [A/F]
gNa = gNa_base * (1 - 0.1 * mode.AF)
    in [mS/uF]
gNa_base = 10 [mS/uF]
    in [mS/uF]
# Parameters
a11 = 3.802 [1/ms] in [1/ms]
a12 = 0.1027
a13 = 2.5 [mV] in [mV]
a14 = 17 [mV] in [mV]
a15 = 0.20
a16 = 150 [mV] in [mV]
a24 = 15 [mV] in [mV]
a25 = 0.23
a34 = 12 [mV] in [mV]
a35 = 0.25
a41 = 0.188495 [ms] in [ms]
a42 = 16.6 [mV] in [mV]
a43 = 0.393956 [ms] in [ms]
a44 = 7 [mV] in [mV]
a51 = 7e-7 [1/ms] in [1/ms]
a52 = 7.2 [mV] in [mV]
a61 = 100
a71 = 0.487e-4 [1/ms] in [1/ms]
a72 = 23.2696 [mV] in [mV]
a81 = 0.1e-7 [1/ms] in [1/ms]
b11 = 0.1917 [1/ms] in [1/ms]
b12 = 20.3 [mV] in [mV]
b21 = 0.2 [1/ms] in [1/ms]
b22 = 2.5 [mV] in [mV]
b31 = 0.22 [1/ms] in [1/ms]
b32 = 7.5 [mV] in [mV]
b51 = 0.0044 [1/ms] in [1/ms]
b52 = 2e-5 [1/ms/mV] in [1/ms/mV]
b61 = 8.9554e-7 [1/ms] in [1/ms]
b62 = 11.3944 [mV] in [mV]
b71 = 0.2868e-3 [1/ms] in [1/ms]
b72 = 35.9898 [mV] in [mV]
b81 = 9.8e-3 [1/ms] in [1/ms]
# Temperature adjustment
Q10_INa = 2.1
Q = Q10_INa^((phys.T - 300 [K]) / 10 [K])
# Rates
a1 = Q * a11 / (a12 * exp(-(V + a13) / a14) + a15 * exp(-(V + a13) / a16))
    in [1/ms]
a2 = Q * a11 / (a12 * exp(-(V + a13) / a24) + a25 * exp(-(V + a13) / a16))
    in [1/ms]
a3 = Q * a11 / (a12 * exp(-(V + a13) / a34) + a35 * exp(-(V + a13) / a16))
    in [1/ms]
a4 = Q / (a41 * exp(-(V + a44) / a42) + a43)
    in [1/ms]
a5 = Q * a51 * exp(-(V + a44) / a52)
    in [1/ms]
a6 = a4 / a61
    in [1/ms]
a7 = Q * a71 * exp(V / a72)
    in [1/ms]
a8 = Q * a81
    in [1/ms]
b1 = Q * b11 * exp(-(V + a13) / b12)
    in [1/ms]
b2 = Q * b21 * exp(-(V - b22) / b12)
    in [1/ms]
b3 = Q * b31 * exp(-(V - b32) / b12)
    in [1/ms]
b4 = (a3 * a4 * a5) / (b3 * b5)
    in [1/ms]
b5 = Q * (b51 + b52 * (V + a44))
    in [1/ms]
b6 = Q * b61 * exp(-V / b62)
    in [1/ms]
b7 = Q * b71 * exp(-V / b72)
    in [1/ms]
b8 = Q * b81
    in [1/ms]
# States
dot(IC3) = b1 * IC2 + b5 * C3 - (a1 + a5) * IC3
dot(IC2) = a1 * IC3 + b2 * IF + b5 * C2 - (b1 + a2 + a5) * IC2
dot(IF) = a4 * O + b5 * C1 + b6 * IM1 + a2 * IC2 - (b4 + a5 + a6 + b2) * IF
dot(IM1) = a6 * IF + b7 * IM2 - (b6 + a7) * IM1
IM2 = 1 - (O + C1 + C2 + C3 + IF + IM1 + IC2 + IC3 + LO + LC1 + LC2 + LC3)
dot(C3) = b8 * LC3 + b1 * C2 + a5 * IC3 - (a1 + b5 + a8) * C3
dot(C2) =  b8 * LC2 + a1 * C3 + b2 * C1 + a5 * IC2 - (b1 + a2 + b5 + a8) * C2
dot(C1) =  b8 * LC1 + a2 * C2 + b3 * O + a5 * IF - (b2 + a3 + b5 + a8) * C1
dot(O) = b8 * LO + a3 * C1 + b4 * IF - (b3 + a4 + a8) * O
dot(LC3) = a8 * C3 + b1 * LC2 - (b8 + a1) * LC3
dot(LC2) = b2 * LC1 + a8 * C2 + a1 * LC3 - (b8 + b1 + a2) * LC2
dot(LC1) = a8 * C1 + a2 * LC2 + b3 * LO - (b8 + b2 + a3) * LC1
dot(LO) = a3 * LC1 + a8 * O - (b8 + b3) * LO

#
# Background sodium current
#
[inab]
use membrane.V
INaB = INaB_jn + INaB_sl
    in [A/F]
INaB_jn = cell.fjn * gNaB * (V - rev.ENa_jn)
    in [A/F]
INaB_sl = cell.fsl * gNaB * (V - rev.ENa_sl)
    in [A/F]
gNaB = 0.597e-3 [mS/uF]
    in [mS/uF]

#
# Na/K pump current
#
[inak]
use membrane.V, phys.FRT
use ion.K_o, sodium.Na_jn, sodium.Na_sl
INaK = INaK_jn + INaK_sl
    in [A/F]
    label I_NaK
INaK_jn = cell.fjn * IbarNaK * fnak * K_o / (1 + (KmNaip / Na_jn)^4) / (K_o + KmKo)
    in [A/F]
INaK_sl = cell.fsl * IbarNaK * fnak * K_o / (1 + (KmNaip / Na_sl)^4) / (K_o + KmKo)
    in [A/F]
IbarNaK = 1.26 [A/F]
    in [A/F]
fnak = 1 / (1 + 0.1245 * exp(-0.1 * V * FRT) + 0.0365 * sigma * exp(-V * FRT))
    sigma = (exp(ion.Na_o / 67.3 [mM]) - 1) / 7
KmKo = 1.5 [mM]
    in [mM]
KmNaip = 11 [mM] * (1 - 0.25 * mode.ISO)
    in [mM]

#
# Rapidly activating potassium current
#
[ikr]
use membrane.V
IKr = gKr * sqrt(ion.K_o / 5.4 [mM]) * x * r * (V - rev.EK)
    in [A/F]
    label I_Kr
gKr = 0.035 [mS/uF]
    in [mS/uF]
    label g_Kr
dot(x) = (inf - x) / tau
    inf = 1 / (1 + exp(-(V + 10 [mV]) / 5 [mV]))
    tau = (+ 550 [ms] / (1 + exp((-22 [mV] - V) / 9 [mV])) * 6 / (1 + exp((V + 11 [mV]) / 9 [mV]))
           + 230 [ms] / (1 + exp((V + 40 [mV]) / 20 [mV])))
        in [ms]
r = 1 / (1 + exp((V + 74 [mV]) / 24 [mV]))

#
# Slowly activating potassium current
#
[iks]
use membrane.V
IKs = IKs_jn + IKs_sl
    in [A/F]
IKs_jn = cell.fjn * gKs * x^2 * (V - rev.EKs)
    in [A/F]
IKs_sl = cell.fsl * gKs * x^2 * (V - rev.EKs)
    in [A/F]
gKs = (1 + mode.AF + 2 * mode.ISO) * gKs_base
    in [mS/uF]
gKs_base = 0.0035 [mS/uF]
    in [mS/uF]
dot(x) = (inf - x) / tau
    inf = 1 / (1 + exp(-(V + 40 [mV] * mode.ISO + 3.8 [mV]) / 14.25 [mV]))
    tau = 990.1 [ms] / (1 + exp(-(V + 40 [mV] * mode.ISO + 2.436 [mV]) / 14.12 [mV]))
        in [ms]

#
# Plateau potassium current
#
[ikp]
use membrane.V
IKp = IKp_jn + IKp_sl
    in [A/F]
IKp_jn = cell.fjn * gKp * kp * (V - rev.EK)
    in [A/F]
IKp_sl = cell.fsl * gKp * kp * (V - rev.EK)
    in [A/F]
gKp = 0.002 [mS/uF]
    in [mS/uF]
kp = 1 / (1 + exp(7.488 - V / 5.98 [mV]))

#
# Muscarinic receptor-activated potassium Current
#
[ikach]
use membrane.V, mode.ACh
IKACh = if(ACh == 0 [mM], 0 [A/F],
           1 / (1 + (0.03 [mM] / ACh)^2.1) * (0.08 [mS/uF] + 0.4 [mS/uF] / (1 + exp((V + 91 [mV]) / 12 [mV]))) * (V - rev.EK))
    in [A/F]

#
# Transient outward potassium current
#
[ito]
use membrane.V
Ito = gto * x * y * (V - rev.EK)
    in [A/F]
    label I_to
gto = (1 - 0.7 * mode.AF) * gto_base
    in [mS/uF]
gto_base = 0.165 [mS/uF]
    in [mS/uF]
dot(x) = (inf - x) / tau
    inf = 1 / (1 + exp(-(V + 1 [mV]) / 11 [mV]))
    tau = 3.5 [ms] * exp(-(V / 30 [mV])^2) + 1.5 [ms]
        in [ms]
dot(y) = (inf - y) / tau
    inf = 1 / (1 + exp((V + 40.5 [mV]) / 11.5 [mV]))
    tau = 25.635 [ms] * exp(-((V + 52.45 [mV]) / 15.8827 [mV])^2) + 24.14 [ms]
        in [ms]

#
# Ultra rapid delayed rectifier outward potassium current
#
[ikur]
use membrane.V
IKur = gKur * O * (V - rev.EK)
    in [A/F]
    label I_Kur
gKur = gKur_base * (1 - 0.5 * mode.AF) * (1 + 2 * mode.ISO)
    in [mS/uF]
gKur_base = 0.045 [mS/uF] * 1.36
    in [mS/uF]
    label g_Kur
a = 2.65 [1/ms] * exp((V - 75.9 [mV]) / 70 [mV])
    in [1/ms]
b = 2.65 [1/ms] * exp((V - 75.9 [mV]) / 70 [mV]) * exp(-(V + 20.9 [mV]) / 10 [mV]) / (1 + 0.3 * exp(-(V + 20.9 [mV]) / 10 [mV]))
    in [1/ms]
kf = 0.0005 [1/ms]
    in [1/ms]
kb = 0.0005 [1/ms] / (1 + exp((V + 10 [mV]) / 10 [mV]))
    in [1/ms]
dot(C1) = b * C2 - 4 * a * C1
dot(C2) = 4 * a * C1 + 2 * b * C3 - C2 * (b + 3 * a)
dot(C3) = 3 * a * C2 + 3 * b * C4 - C3 * (2 * b + 2 * a)
dot(C4) = 2 * a * C3 + 4 * b * O - C4 * (3 * b + a)
dot(O) = a * C4 + kb * I - O * (4 * b + kf)
I = 1 - C1 - C2 - C3 - C4 - O

#
# Inward recitifer potassium current
#
[ik1]
use membrane.V, rev.EK
IK1 = gK1 * inf * (V - EK)
    in [A/F]
    label I_K1
gK1 = (1 + mode.AF) * sqrt(ion.K_o / 5.4 [mM]) * gK1_base
    in [mS/uF]
gK1_base = 0.0525 [mS/uF]
    in [mS/uF]
inf = a / (a + b)
    a = 1.02 / (1 + exp(0.2385 [1/mV] * (V - EK - 59.215 [mV])))
    b = (0.49124 * exp(0.08032 [1/mV] * (V - EK + 5.476 [mV]))
         + exp(0.06175 [1/mV] * (V - EK - 594.31 [mV]))
        ) / (1 + exp(-0.5143 [1/mV] * (V - EK + 4.753 [mV])))

#
# Calcium activated Chloride current
#
[iclca]
use membrane.V
IClCa = IClCa_jn + IClCa_sl
    in [A/F]
IClCa_jn = cell.fjn * gClCa / (1 + KdClCa / calcium.Ca_jn) * (V - rev.ECl)
    in [A/F]
IClCa_sl = cell.fsl * gClCa / (1 + KdClCa / calcium.Ca_sl) * (V - rev.ECl)
    in [A/F]
gClCa = 0.0548 [mS/uF]
    in [mS/uF]
KdClCa = 0.1 [mM]
    in [mM]

#
# Background chloride current
#
[iclb]
IClB = gClB * (membrane.V - rev.ECl)
    in [A/F]
gClB = 9e-3 [mS/uF]
    in [mS/uF]

#
# L-type calcium current
#
[ical]
use membrane.V
use phys.F, phys.FRT
use calcium.Ca_jn, calcium.Ca_sl, ion.Ca_o
use sodium.Na_jn, sodium.Na_sl, ion.Na_o
use ion.K_i, ion.K_o
use cell.fjn_CaL, cell.fsl_CaL
use mode.ISO, mode.AF
ICaL = ICaL_Ca + ICaL_Na + ICaL_K
    in [A/F]
    label I_CaL
ICaL_Ca = ICaL_Ca_jn + ICaL_Ca_sl
    in [A/F]
ICaL_Ca_jn = fjn_CaL * ibarca_jn * d * f * (1 - fCa_jn) * Q
    in [A/F]
ICaL_Ca_sl = fsl_CaL * ibarca_sl * d * f * (1 - fCa_sl) * Q
    in [A/F]
ICaL_Na = ICaL_Na_jn + ICaL_Na_sl
    in [A/F]
ICaL_Na_jn = fjn_CaL * ibarna_jn * d * f * (1 - fCa_jn) * Q
    in [A/F]
ICaL_Na_sl = fsl_CaL * ibarna_sl * d * f * (1 - fCa_sl) * Q
    in [A/F]
ICaL_K = ibark * d * f * (fjn_CaL * (1 - fCa_jn) + cell.fsl_CaL * (1 - fCa_sl)) * Q
    in [A/F]
# Permeability
scaling = 1
    label g_CaL
PCa = scaling * (1 + 0.5 * ISO) * (1 - 0.5 * AF) * 2.7e-4 [L/F/ms]
    in [L/F/ms]
PNa = scaling * (1 + 0.5 * ISO) * (1 - 0.5 * AF) * 0.75e-8 [L/F/ms]
    in [L/F/ms]
PK = scaling * (1 + 0.5 * ISO) * (1 - 0.5 * AF) * 1.35e-7 [L/F/ms]
    in [L/F/ms]
Q10CaL = 1.8
Q = Q10CaL^phys.Q * 0.45
# Driving terms
ibarca_jn = PCa * 4 * V * F * FRT * (0.341 * Ca_jn * exp(2 * V * FRT) - 0.341 * Ca_o) / (exp(2 * V * FRT) - 1)
    in [A/F]
ibarca_sl = PCa * 4 * V * F * FRT * (0.341 * Ca_sl * exp(2 * V * FRT) - 0.341 * Ca_o) / (exp(2 * V * FRT) - 1)
    in [A/F]
ibarna_jn = PNa * V * F * FRT * (0.75 * Na_jn * exp(V * FRT) - 0.75 * Na_o) / (exp(V * FRT) - 1)
    in [A/F]
ibarna_sl = PNa * V * F * FRT * (0.75 * Na_sl * exp(V * FRT) - 0.75 * Na_o) / (exp(V * FRT) - 1)
    in [A/F]
ibark = PK * V * F * FRT * (0.75 * K_i * exp(V * FRT) - 0.75 * K_o) / (exp(V * FRT) - 1)
    in [A/F]
# Activation
dot(d) = (inf - d) / tau
    inf = 1 / (1 + exp(-(V + 3 [mV] * ISO + 9 [mV]) / 6 [mV]))
    tau = inf * (1 - exp(-(V + 3 [mV] * ISO + 9 [mV]) / 6 [mV])) / (0.035 [1/ms/mV] * (V + 3 [mV] * ISO + 9 [mV]))
        in [ms]
# Voltage-dependent inactivation
dot(f) = (inf - f) / tau
    inf = 1 / (1 + exp((V + 3 [mV] * ISO + 30 [mV]) / 7 [mV])) + 0.2 / (1 + exp((50 [mV] - V - 3 [mV] * ISO) / 20 [mV]))
    tau = 1 / (0.0197 [1/ms] * exp(-(0.0337 [1/mV] * (V + 3 [mV] * ISO + 25 [mV]))^2) + 0.02 [1/ms])
        in [ms]
# Calcium-dependent inactivation (gate = 1 minus these)
dot(fCa_jn) = 1.7 [1/mM/ms] * Ca_jn * (1 - fCa_jn) - 11.9e-3 [1/ms] * fCa_jn
dot(fCa_sl) = 1.7 [1/mM/ms] * Ca_sl * (1 - fCa_sl) - 11.9e-3 [1/ms] * fCa_sl

#
# Sodium/Calcium exchanger current
#
[inaca]
use membrane.V
use phys.FRT
use sodium.Na_jn, sodium.Na_sl, ion.Na_o
use calcium.Ca_jn, calcium.Ca_sl, ion.Ca_o
INaCa = INaCa_jn + INaCa_sl
    in [A/F]
    label I_NaCa
INaCa_jn = cell.fjn * IbarNaCa * Q * Ka * (s1 - s2) / s3 / (1 + ksat * exp((nu - 1) * V * FRT))
    in [A/F]
    s1 = exp(nu * V * FRT) * Na_jn^3 * Ca_o
        in [mM^4]
    s2 = exp((nu - 1) * V * FRT) * Na_o^3 * Ca_jn
        in [mM^4]
    s3 = KmCai * Na_o^3 * (1 + (Na_jn / KmNai)^3) + KmNao^3 * Ca_jn * (1 + Ca_jn / KmCai) + KmCao * Na_jn^3 + Na_jn^3 * Ca_o + Na_o^3 * Ca_jn
        in [mM^4]
    Ka = 1 / (1 + (Kdact / Ca_jn)^2)
INaCa_sl = cell.fsl * IbarNaCa * Q * Ka * (s1 - s2) / s3 / (1 + ksat * exp((nu - 1) * V * FRT))
    in [A/F]
    s1 = exp(nu * V * FRT) * Na_sl^3 * Ca_o
        in [mM^4]
    s2 = exp((nu - 1) * V * FRT) * Na_o^3 * Ca_sl
        in [mM^4]
    s3 = KmCai * Na_o^3 * (1 + (Na_sl / KmNai)^3) + KmNao^3 * Ca_sl * (1 + Ca_sl / KmCai) + KmCao * Na_sl^3 + Na_sl^3 * Ca_o + Na_o^3 * Ca_sl
        in [mM^4]
    Ka = 1 / (1 + (Kdact / Ca_sl)^2)
IbarNaCa = (1 + 0.4 * mode.AF) * 3.15 [A/F]
    in [A/F]
Q10NaCa = 1.57
Q = Q10NaCa^phys.Q
Kdact = 0.384e-3 [mM]
    in [mM]
KmNai = 12.29 [mM]
    in [mM]
KmNao = 87.5 [mM]
    in [mM]
KmCai = 3.59e-3 [mM]
    in [mM]
KmCao = 1.3 [mM]
    in [mM]
ksat = 0.27
nu = 0.35

#
# Sarcolemmal calcium pump current
#
[ipca]
use calcium.Ca_jn, calcium.Ca_sl
IpCa = IpCa_jn + IpCa_sl
    in [A/F]
IpCa_jn = cell.fjn * Q * IbarSLCaP * a / (a + b)
    in [A/F]
    a = (calcium.Ca_jn / 1 [mM])^1.6
IpCa_sl = cell.fsl * Q * IbarSLCaP * a / (a + b)
    in [A/F]
    a = (calcium.Ca_sl / 1 [mM])^1.6
b = (KmPCa / 1 [mM])^1.6
KmPCa = 0.0005 [mM]
    in [mM]
IbarSLCaP = 0.0471 [A/F]
    in [A/F]
Q10 = 2.35
Q = Q10^phys.Q

#
# Background calcium current
#
[icab]
use membrane.V
ICaB = ICaB_jn + ICaB_sl
    in [A/F]
ICaB_jn = cell.fjn * gCaB * (V - rev.ECa_jn)
    in [A/F]
ICaB_sl = cell.fsl * gCaB * (V - rev.ECa_sl)
    in [A/F]
gCaB = 6.0643e-4 [mS/uF]
    in [mS/uF]

#
# SR calcium release
#
[ryr]
use calcium.Ca_jn, cabsr.Ca_sr
J_SRCarel = ks * o * (Ca_sr - Ca_jn)
    in [mM/ms]
ks = 25 [1/ms]
    in [1/ms]
kom = 0.06 [1/ms]
    in [1/ms]
kim = 0.005 [1/ms]
    in [1/ms]
kiCa = 0.5 [1/mM/ms]
    in [1/mM/ms]
koCa = 10 [1/mM^2/ms] * (1 + 2 * mode.AF + mode.ISO * (1 - mode.AF))
    in [1/mM^2/ms]
ec50SR = 0.45 [mM]
    in [mM]
MaxSR = 15
MinSR = 1
kCaSR = MaxSR - (MaxSR - MinSR) / (1 + (ec50SR / Ca_sr)^2.5)
koSRCa = koCa / kCaSR
    in [1/mM^2/ms]
kiSRCa = kiCa * kCaSR
    in [1/mM/ms]
dot(r) = (kim * ri - kiSRCa * Ca_jn * r) - (koSRCa * Ca_jn^2 * r - kom * o)
dot(o) = (koSRCa * Ca_jn^2 * r - kom * o) - (kiSRCa * Ca_jn * o - kim * i)
dot(i) = (kiSRCa * Ca_jn * o - kim * i) - (kom * i - koSRCa * Ca_jn^2 * ri)
ri = 1 - r - o - i

#
# SR calcium pump
#
[serca]
use calcium.Ca_i, cabsr.Ca_sr
J_serca = Q * Vmax * ((Ca_i / Kmf)^h - (Ca_sr / Kmr)^h) / (1 + (Ca_i / Kmf)^h + (Ca_sr / Kmr)^h)
    in [mM/ms]
Vmax = 5.3114e-3 [mM/ms]
    in [mM/ms]
Kmr = 1.7 [mM]
    in [mM]
Kmf = 0.246e-3 [mM] * (2.5 - 1.5 * mode.ISO)  # 1.5 was 1.25 in Grandi 2011
    in [mM]
h = 1.787
Q = Q10^phys.Q
Q10 = 2.6

#
# SR calcium leak
#
[caleak]
use cabsr.Ca_sr, calcium.Ca_jn
J_SRleak = 5.348e-6 [1/ms] * (Ca_sr - Ca_jn) * (1 + 0.25 * mode.AF)
    in [mM/ms]

#
# Cytosolic calcium buffering
#
[cabct]
use cell.Vmyo, cell.Vjn, cell.Vsl
use calcium.Ca_i, calcium.Ca_jn, calcium.Ca_sl
use cabss.J_CaB_jn, cabss.J_CaB_sl
use ion.Mg_i
# In Grandi 2011, this had additional terms dot(TnCHMg) + dot(MyoMg)
J_CaB_cytosol = dot(TnCL) + dot(TnCHCa) + dot(CaM) + dot(MyoCa) + dot(SRB)
    in [mM/ms]
# TnCl
Bmax_TnClow = 0.07 [mM]
    in [mM]
koff_tncl = 19.6e-3 [1/ms] * (1 + 0.5 * mode.ISO)
    in [1/ms]
kon_tncl = 32.7 [1/mM/ms]
    in [1/mM/ms]
Bmax_TnChigh = 140e-3 [mM]
    in [mM]
koff_tnchCa = 0.032e-3 [1/ms]
    in [1/ms]
kon_tnchCa = 2.37 [1/mM/ms]
    in [1/mM/ms]
koff_tnchMg = 3.33e-3 [1/ms]
    in [1/ms]
kon_tnchMg = 3e-3 [1/mM/ms]
    in [1/mM/ms]
dot(TnCL) = kon_tncl * Ca_i * (Bmax_TnClow - TnCL) - koff_tncl * TnCL
    in [mM]
dot(TnCHCa) = kon_tnchCa * Ca_i * (Bmax_TnChigh - TnCHCa - TnCHMg) - koff_tnchCa * TnCHCa
    in [mM]
dot(TnCHMg) = kon_tnchMg * Mg_i * (Bmax_TnChigh - TnCHCa - TnCHMg) - koff_tnchMg * TnCHMg
    in [mM]
# Calmodulin
Bmax_CaM = 24e-3 [mM]
    in [mM]
koff_cam = 238e-3 [1/ms]
    in [1/ms]
kon_cam = 34 [1/mM/ms]
    in [1/mM/ms]
dot(CaM) = kon_cam * Ca_i * (Bmax_CaM - CaM) - koff_cam * CaM
    in [mM]
# Myosin-Ca and Myosin-Mg
Bmax_myosin = 140e-3 [mM]
    in [mM]
koff_myoCa = 0.46e-3 [1/ms]
    in [1/ms]
kon_myoCa = 13.8 [1/mM/ms]
    in [1/mM/ms]
koff_myoMg = 0.057e-3 [1/ms]
    in [1/ms]
kon_myoMg = 0.0157 [1/mM/ms]
    in [1/mM/ms]
dot(MyoCa) = kon_myoCa * Ca_i * (Bmax_myosin - MyoCa - MyoMg) - koff_myoCa * MyoCa
    in [mM]
dot(MyoMg) = kon_myoMg * Mg_i * (Bmax_myosin - MyoCa - MyoMg) - koff_myoMg * MyoMg
    in [mM]
# SRB
Bmax_SR = 19 * 0.9e-3 [mM]
    in [mM]
koff_sr = 60e-3 [1/ms]
    in [1/ms]
kon_sr = 100 [1/mM/ms]
    in [1/mM/ms]
dot(SRB) = kon_sr * Ca_i * (Bmax_SR - SRB) - koff_sr * SRB
    in [mM]

#
# Junctional and SL calcium buffers
#
[cabss]
use calcium.Ca_sl, calcium.Ca_jn
use cell.Vmyo, cell.Vsl, cell.Vjn
J_CaB_jn = dot(SLL_jn) + dot(SLH_jn)
    in [mM/ms]
J_CaB_sl = dot(SLL_sl) + dot(SLH_sl)
    in [mM/ms]
Bmax_SLlowjn = 4.6e-3 [mM] * Vmyo / Vjn * 0.1
    in [mM]
Bmax_SLlowsl = 37.4e-3 [mM] * Vmyo / Vsl
    in [mM]
Bmax_SLhighjn = 1.65e-3 [mM] * Vmyo / Vjn * 0.1
    in [mM]
Bmax_SLhighsl = 13.4e-3 [mM] * Vmyo / Vsl
    in [mM]
koff_sll = 1.3 [1/ms]
    in [1/ms]
kon_sll = 100 [1/mM/ms]
    in [1/mM/ms]
koff_slh = 0.03 [1/ms]
    in [1/ms]
kon_slh = 100 [1/mM/ms]
    in [1/mM/ms]
dot(SLL_jn) = kon_sll * Ca_jn * (Bmax_SLlowjn - SLL_jn) - koff_sll * SLL_jn
    in [mM]
dot(SLL_sl) = kon_sll * Ca_sl * (Bmax_SLlowsl - SLL_sl) - koff_sll * SLL_sl
    in [mM]
dot(SLH_jn) = kon_slh * Ca_jn * (Bmax_SLhighjn - SLH_jn) - koff_slh * SLH_jn
    in [mM]
dot(SLH_sl) = kon_slh * Ca_sl * (Bmax_SLhighsl - SLH_sl) - koff_slh * SLH_sl
    in [mM]

#
# SR calcium buffer
#
[cabsr]
use cell.Vmyo, cell.Vsr
koff_csqn = 65 [1/ms]
    in [1/ms]
kon_csqn = 100 [1/mM/ms]
    in [1/mM/ms]
Bmax_Csqn = 0.14 [mM] * Vmyo / Vsr
    in [mM]
dot(Csqn) = kon_csqn * Ca_sr * (Bmax_Csqn - Csqn) - koff_csqn * Csqn
    in [mM]
dot(Ca_sr) = serca.J_serca - (caleak.J_SRleak * Vmyo / Vsr + ryr.J_SRCarel) - dot(Csqn)
    in [mM]

#
# Calcium concentrations
#
[calcium]
use cell.C, phys.F
use cell.Vjn, cell.Vsl, cell.Vmyo, cell.Vsr
use serca.J_serca, ryr.J_SRCarel, caleak.J_SRleak
use cabct.J_CaB_cytosol, cabss.J_CaB_jn, cabss.J_CaB_sl
J_Ca_jnsl = 1 [L/ms] / 1.2134e12
    in [L/ms]
J_Ca_slmyo = 1 [L/ms] / 2.68510e11
    in [L/ms]
I_Ca_tot_jn = ical.ICaL_Ca_jn + icab.ICaB_jn + ipca.IpCa_jn - 2 * inaca.INaCa_jn
    in [A/F]
I_Ca_tot_sl = ical.ICaL_Ca_sl + icab.ICaB_sl + ipca.IpCa_sl - 2 * inaca.INaCa_sl
    in [A/F]
dot(Ca_i) = -J_serca * Vsr / Vmyo - J_CaB_cytosol + J_Ca_slmyo / Vmyo * (Ca_sl - Ca_i)
    in [mM]
    label Ca_i
dot(Ca_jn) = -I_Ca_tot_jn * C / (Vjn * 2 * F) + J_Ca_jnsl / Vjn * (Ca_sl - Ca_jn) - J_CaB_jn + J_SRCarel * Vsr / Vjn + J_SRleak * Vmyo / Vjn
    in [mM]
dot(Ca_sl) = -I_Ca_tot_sl * C / (Vsl * 2 * F) + J_Ca_jnsl / Vsl * (Ca_jn - Ca_sl) + J_Ca_slmyo / Vsl * (Ca_i - Ca_sl) - J_CaB_sl
    in [mM]

#
# Sodium buffering
#
[nabss]
Bmax_Na_jn = 7.561 [mM]
    in [mM]
Bmax_Na_sl = 1.65 [mM]
    in [mM]
koff_na = 1e-3 [1/ms]
    in [1/ms]
kon_na = 1e-4 [1/mM/ms]
    in [1/mM/ms]
dot(NaB_jn) = kon_na * sodium.Na_jn * (Bmax_Na_jn - NaB_jn) - koff_na * NaB_jn
    in [mM]
dot(NaB_sl) = kon_na * sodium.Na_sl * (Bmax_Na_sl - NaB_sl) - koff_na * NaB_sl
    in [mM]

#
# Sodium concentrations
#
[sodium]
use cell.C, phys.F
use cell.Vjn, cell.Vsl, cell.Vmyo
J_Na_jnsl = 1 [L/ms] / (1.6382e12 / 3 * 100)
    in [L/ms]
J_Na_slmyo = 1 [L/ms] / (1.8308e10 / 3 * 100)
    in [L/ms]
I_Na_tot_jn = ina.INa_jn + inab.INaB_jn + 3 * inaca.INaCa_jn + 3 * inak.INaK_jn + ical.ICaL_Na_jn
    in [A/F]
I_Na_tot_sl = ina.INa_sl + inab.INaB_sl + 3 * inaca.INaCa_sl + 3 * inak.INaK_sl + ical.ICaL_Na_sl
    in [A/F]
dot(Na_jn) = -I_Na_tot_jn * C / (Vjn * F) + J_Na_jnsl / Vjn * (Na_sl - Na_jn) - dot(nabss.NaB_jn)
    in [mM]
dot(Na_sl) = -I_Na_tot_sl * C / (Vsl * F) + J_Na_jnsl / Vsl * (Na_jn - Na_sl) + J_Na_slmyo / Vsl * (Na_i - Na_sl) - dot(nabss.NaB_sl)
    in [mM]
dot(Na_i) = J_Na_slmyo / Vmyo * (Na_sl - Na_i)
    in [mM]
    label Na_i

#
# Potassium concentration
#
[potassium]
I_K_tot = ito.Ito + ikr.IKr + iks.IKs + ik1.IK1 - 2 * inak.INaK + ical.ICaL_K + ikp.IKp + ikur.IKur + ikach.IKACh
    in [A/F]

[[protocol]]
# Level  Start    Length   Period   Multiplier
1        50       0.5      1000     0

[[script]]
import matplotlib.pyplot as plt
import myokit

# Get model and protocol, create simulation
m = get_model()
p = get_protocol()
s = myokit.Simulation(m, p)

# Run simulation
s.pre(10*1000)
d = s.run(2000, log=['engine.time', 'membrane.V', 'calcium.Ca_i'])

d = d.npview()
d['engine.time'] /= 1e3

# Display the results
plt.figure(figsize=(9, 7))
ax = plt.subplot(2, 1, 1)
ax.plot(d.time(), d['membrane.V'])
ax = plt.subplot(2, 1, 2)
ax.plot(d.time(), d['calcium.Ca_i'], label='Ca_i')
plt.show()