fabbri-2017.mmt

[[model]]
name: fabbri-2017
version: 20240904
mmt_authors: Michael Clerx, Aditi Agrawal
display_name: Fabbri et al., 2017
desc: """
    Model of human sino-atrial node (SAN) cells by Fabbri et al. [1], based on
    an earlier rabbit model by Severi et al. [2].

    This Myokit implementation is based on the CellML code by Alan Fabbri [3],
    which was tested in a reproducibility study [4]. The Myokit implementation
    has extensive reformatting and unit corrections, and includes a small
    offset in IKACh to avoid divide-by-zero errors, as suggested by Alan
    Fabbri. It was tested against the original CellML by comparing the
    calculated derivatives, which matched to within machine precision.

    [1] Fabbri, A. Fantini, M., Wilders, R., & Severi, S. (2017) Computational
        analysis of the human sinus node action potential: model development
        and effects of mutations. Journal of Physiology, 595.7 pp 2365-2396.
        https://doi.org/10.1113/JP273259

    [2] Severi, S., Fantini, M., Charawi, L. A., & DiFrancesco, D. (2012) An
        updated computational model of rabbit sinoatrial action potential to
        investigate the mechanisms of heart rate modulation. Journal of
        Physiology 590.18 pp 4483-4499.
        https://doi.org/10.1113/jphysiol.2012.229435

    [3] CellML code by Alan Fabbri
        https://www.mcbeng.it/en/downloads/software/hap-san.html

    [4] Physiome reproduction paper
        Afshar, N., Fabbri, A., Severi, S., Garny, A., Nickerson, D. (2021)
        Reproducibility study of the Fabbri et al. 2017 model of the human
        sinus node action potential. Physiome.
        https://doi.org/10.36903/physiome.16550526

"""
# Initial values
membrane.V     = -4.7787168e1
ina.m          = 0.447724
ina.h          = 0.003058
ical.d         = 0.001921
ical.f         = 0.846702
ical.fCa       = 0.844449
icat.d         = 0.268909
icat.f         = 0.020484
if.y           = 0.009508
ikur.r         = 0.011845
ikur.s         = 0.845304
ito.q          = 0.430836
ito.r          = 0.014523
ikr.pas        = 0.283185
ikr.paf        = 0.011068
ikr.piy        = 0.709051
iks.n          = 0.1162
ikach.a        = 0.00277
carel.r        = 0.9308
carel.o        = 6.181512e-9
carel.i        = 4.595622e-10
carel.ri       = 0.069199
sodium.Na_i    = 5
calcium.Ca_i    = 9.15641e-6
calcium.Ca_sub = 6.226104e-5
calcium.Ca_nsr = 0.435148
calcium.Ca_jsr = 0.409551
cabuf.fTC      = 0.017929
cabuf.fTMC     = 0.259947
cabuf.fTMM     = 0.653777
cabuf.fCMi     = 0.217311
cabuf.fCMs     = 0.158521
cabuf.fCQ      = 0.138975

#
# Simulation engine
#
[engine]
time = 0 [s] bind time
    in [s]

#
# Membrane potential
#
[membrane]
dot(V) = -i_tot / cell.C
    in [mV]
    label membrane_potential
i_tot = (
        + if.If
        + ikr.IKr
        + iks.IKs
        + ito.Ito
        + inak.INaK
        + inaca.INaCa
        + ina.INa
        + ical.ICaL
        + icat.ICaT
        + ikach.IKACh
        + ikur.IKur
    )
    in [nA]
    label cellular_current

#
# Acetylcholine
#
[ach]
ACh = 0 [mM]
    in [mM]
    desc: Acetylcholine concentration

#
# Isoproterenol
#
[iso]
iso = 0
    desc: Set to 1 to simulate with 1 uM of isoproterenol present

#
# Cell geometry
#
[cell]
C = 5.7e-5 [uF]
    in [uF]
    label membrane_capacitance
L = 67 [um]
    in [um]
R = 3.9 [um]
    in [um]
pi = 3.14159265358979312
vcell = 1e-9 [uL/um^3] * pi * R^2 * L
    in [uL]
Lsub = 0.02 [um]
    in [um]
vsub = 1e-9 [uL/um^3] * 2 * pi * Lsub * (R - Lsub / 2) * L
    in [uL]
vi = 0.46 * vcell - vsub
    in [uL]
vjsr = 0.0012 * vcell
    in [uL]
vnsr = 0.0116 * vcell
    in [uL]

#
# Physical constants and temperature
#
[phys]
F =  9.64853415e4 [C/mol]
    in [C/mol]
R = 8314.472 [J/kmol/K]
    in [mJ/mol/K]
RTF = R * T / F
    in [mV]
T = 310 [K]
    in [K]

#
# Extracellular concentrations
#
[extra]
Ca_o = 1.8 [mM]
    in [mM]
Na_o = 140 [mM]
    in [mM]
K_o = 5.4 [mM]
    in [mM]

#
# Reversal potentials
#
[rev]
E_K = phys.RTF * log(extra.K_o / potassium.K_i)
    in [mV]
E_Na = phys.RTF * log(extra.Na_o / sodium.Na_i)
    in [mV]
PKNa = 0.12
E_Ks = phys.RTF * log((extra.K_o + PKNa * extra.Na_o) / (potassium.K_i + PKNa * sodium.Na_i))
    in [mV]
E_mh = phys.RTF * log((extra.Na_o + PKNa * extra.K_o) / (sodium.Na_i + PKNa * potassium.K_i))
    in [mV]

#
# Fast sodium current
#
[ina]
use membrane.V
gNa = 0.0223 [uS]
    in [uS]
gNaL = 0 [uS]
    in [uS]
INaF = gNa * m^3 * h * (V - rev.E_mh)
    in [nA]
INaL = gNaL * m^3 * (V - rev.E_mh)
    in [nA]
INa = INaF + INaL
    in [nA]
dot(m) = (inf - m) / tau
    inf = 1 / (1 + exp(-(V + 42.0504 [mV]) / 8.3106 [mV]))
    tau = 1 / (alpha + beta)
        in [s]
    E0 = V + 41 [mV]
        in [mV]
    alpha = if(abs(E0) < 1e-5 [mV], 2000 [1/s], 200 [1/mV/s] * E0 / (1 - exp(-0.1 [1/mV] * E0)))
        in [1/s]
    beta = 8000 [1/s] * exp(-0.056 [1/mV] * (V + 66 [mV]))
        in [1/s]
dot(h) = (inf - h) / tau
    inf = 1 / (1 + exp((V + 69.804 [mV]) / 4.4565 [mV]))
    tau = 1 / (alpha + beta)
        in [s]
    alpha = 20 [1/s] * exp(-0.125 [1/mV] * (V + 75 [mV]))
        in [1/s]
    beta = 2000 [1/s] / (320 * exp(-0.1 [1/mV] * (V + 75 [mV])) + 1)
        in [1/s]

#
# L-type calcium current
#
[ical]
use membrane.V, phys.RTF
use calcium.Ca_sub, sodium.Na_i, potassium.K_i
use extra.Ca_o, extra.Na_o, extra.K_o
dot(d) = (inf - d) / tau
    inf = 1 / (1 + exp(-(V + 16.4508 [mV] - iso_shift) / (4.3371 [mV] * (1 + iso_slope / 100))))
        iso_shift = piecewise(iso.iso > 0, -8 [mV], 0 [mV])
            in [mV]
        iso_slope = piecewise(iso.iso > 0, -27, 0)
    tau = 0.001 [s/ms] / (alpha + beta)
        in [s]
        av = piecewise(V == -41.8 [mV], -41.80001 [mV], V == -6.8 [mV], -6.80001 [mV], V)
            in [mV]
        bv = piecewise(V == -1.8 [mV], -1.80001 [mV], V)
            in [mV]
        alpha = -0.02839 [1/ms/mV] * (av + 41.8 [mV]) / (exp(-(av + 41.8 [mV]) / 2.5 [mV]) - 1) - 0.0849 [1/mV/ms] * (av + 6.8 [mV]) / (exp(-(av + 6.8 [mV]) / 4.8 [mV]) - 1)
            in [1/ms]
        beta = 0.01143 [1/ms/mV] * (bv + 1.8 [mV]) / (exp((bv + 1.8 [mV]) / 2.5 [mV]) - 1)
            in [1/ms]
dot(f) = (inf - f) / tau
    tau = 0.001 [s/ms] * (44.3 [ms] + 230 [ms] * exp(-((V + 36 [mV]) / 10 [mV])^2))
        in [s]
    inf = 1 / (1 + exp((V + 37.4 [mV]) / 5.3 [mV]))
dot(fCa) = (inf - fCa) / tau
    inf = Km / (Km + calcium.Ca_sub)
    tau = 0.001 [s/ms] * inf / alpha
        in [s]
    Km = 0.000338 [mM]
        in [mM]
    alpha = 0.0075 [1/ms]
        in [1/ms]
# Block and increase
ACh_block = 0.31 * ach.ACh / (ach.ACh + 9e-5 [mM])
Iso_increase = if(iso.iso > 0, 1.23, 1)
# Current
P_CaL = 0.4578 [nA/mM]
    in [nA/mM]
IsiCa = d * f * fCa * V / (RTF * (1 - exp(-V * 2 / RTF))) * (Ca_sub - Ca_o * exp(-2 * V / RTF)) * P_CaL * 2
    in [nA]
IsiK  = d * f * fCa * V / (RTF * (1 - exp(-V / RTF))) * (K_i - K_o * exp(-V / RTF)) * P_CaL * 0.000365
    in [nA]
IsiNa = d * f * fCa * V / (RTF * (1 - exp(-V / RTF))) * (Na_i - Na_o * exp(-V / RTF)) * P_CaL * 1.85e-5
    in [nA]
ICaL = (IsiCa + IsiK + IsiNa) * (1 - ACh_block) * 1 * Iso_increase
    in [nA]

#
# T-type calcium current
#
[icat]
use membrane.V, phys.RTF
P_CaT = 0.04132 [nA/mM]
    in [nA/mM]
ICaT = 2 * P_CaT * V / (RTF * (1 - exp(-1 * V * 2 / RTF))) * (calcium.Ca_sub - extra.Ca_o * exp(-2 * V / RTF)) * d * f
    in [nA]
dot(d) = (inf - d) / tau
    inf = 1 / (1 + exp(-(V + 38.3 [mV]) / 5.5 [mV]))
    tau = 0.001 [s] / (1.068 * exp((V + 38.3 [mV]) / 30 [mV]) + 1.068 * exp(-(V + 38.3 [mV]) / 30 [mV]))
        in [s]
dot(f) = (inf - f) / tau
    inf = 1 / (1 + exp((V + 58.7 [mV]) / 3.8 [mV]))
    tau = 1 [s] / (16.67 * exp(-(V + 75 [mV]) / 83.3 [mV]) + 16.67 * exp((V + 75 [mV]) / 15.38 [mV]))
        in [s]

#
# Transient outward potassium current
#
[ito]
use membrane.V
gto = 0.0035 [uS]
    in [uS]
Ito = gto * (V - rev.E_K) * q * r
    in [nA]
dot(q) = (inf - q) / tau
    inf = 1 / (1 + exp((V + 49 [mV]) / 13 [mV]))
    tau = 0.001 [s/ms] * 0.6 * (65.17 [ms] / (0.57 * exp(-0.08 [1/mV] * (V + 44 [mV])) + 0.065 * exp(0.1 [1/mV] * (V + 45.93 [mV]))) + 10.1 [ms])
        in [s]
dot(r) = (inf - r) / tau
    inf = 1 / (1 + exp(-(V - 19.3 [mV]) / 15 [mV]))
    tau = 0.001 [s/ms] * 0.66 * 1.4 * (15.59 [ms] / (1.037 * exp(0.09 [1/mV] * (V + 30.61 [mV])) + 0.369 * exp(-0.12 [1/mV] * (V + 23.84 [mV]))) + 2.98 [ms])
        in [s]

#
# Rapid delayed-rectifier potassium current
#
[ikr]
use membrane.V
gKr = 0.00424 [uS]
    in [uS]
    label g_Kr
IKr = gKr * (0.9 * paf + 0.1 * pas) * piy * (V - rev.E_K)
    in [nA]
a_inf = 1 / (1 + exp(-(V + 10.0144 [mV]) / 7.6607 [mV]))
dot(paf) = (a_inf - paf) / tau
    tau = 1 [s] / (30 * exp(V / 10 [mV]) + exp(-V / 12 [mV]))
        in [s]
dot(pas) = (a_inf - pas) / tau
    tau =  8.4655354e-1 [s] / (4.2 * exp(V / 17 [mV]) + 0.15 * exp(-V / 21.6 [mV]))
        in [s]
dot(piy) = (inf - piy) / tau
    inf = 1 / (1 + exp((V + 28.6 [mV]) / 17.1 [mV]))
    tau = 1 [s] / (100 * exp(-V / 54.645 [mV]) + 656 * exp(V / 106.157 [mV]))
        in [s]

#
# Slow delayed-rectifier potassium current
#
[iks]
use membrane.V
gKs = 0.00065 [uS] * if(iso.iso > 0, 1.2, 1)
    in [uS]
IKs = gKs * (V - rev.E_Ks) * n^2
    in [nA]
iso_shift = if(iso.iso > 0, -14 [mV], 0 [mV])
    in [mV]
dot(n) = (inf - n) / tau
    inf = sqrt(1 / (1 + exp(-(V + 0.6383 [mV] - iso_shift) / 10.7071 [mV])))
    tau = 1 / (alpha + beta)
        in [s]
    alpha = 28 [1/s] / (1 + exp(-(V - 40 [mV] - iso_shift) / 3 [mV]))
        in [1/s]
    beta = 1 [1/s] * exp(-(V - iso_shift - 5 [mV]) / 25 [mV])
        in [1/s]

#
# Ultra-rapid potassium current
#
[ikur]
use membrane.V
gKur = 0.0001539 [uS]
    in [uS]
IKur = gKur * r * s * (V - rev.E_K)
    in [nA]
dot(r) = (inf - r) / tau
    inf = 1 / (1 + exp((V + 6 [mV]) / -8.6 [mV]))
    tau = 0.009 [s] / (1 + exp((V + 5 [mV]) / 12 [mV])) + 0.0005 [s]
        in [s]
dot(s) = (inf - s) / tau
    inf = 1 / (1 + exp((V + 7.5 [mV]) / 10 [mV]))
    tau = 0.59 [s] / (1 + exp((V + 60 [mV]) / 10 [mV])) + 3.05 [s]
        in [s]

#
# Acetylcholine-sensitive potassium current
#
[ikach]
use membrane.V
use ach.ACh
ACh_on = 1
gKACh = 0.00345 [uS]
    in [uS]
IKACh = if(ACh > 0 [mM],
        a * ACh_on * gKACh * (V - rev.E_K) * (1 + exp((V + 20 [mV]) / 20 [mV])),
        0 [nA])
    in [nA]
dot(a) = alpha * (1 - a) - beta * a
    alpha = 0.025641 [1/s] + alpha_ach
        in [1/s]
        desc: Alpha gate with slight update to avoid divide-by-zero as suggested by Alan Fabbri
    alpha_ach = if(ACh == 0 [mM], 0 [1/s],
                   (3.5988 [1/s] - 0.025641 [1/s]) / (1 + 1.2155e-6 / (ACh / 1 [mM])^1.6951))
        in [1/s]
    beta = 10 [1/s] * exp(0.0133 [1/mV] * (V + 40 [mV]))
        in [1/s]

#
# Funny current
#
[if]
use membrane.V
use extra.K_o
If = IfNa + IfK
    in [nA]
IfK = y * gfK * (V - rev.E_K)
    in [nA]
IfNa = y * gfNa * (V - rev.E_Na)
    in [nA]
gfK = gf / (alpha + 1)
    in [uS]
gfNa = gf * alpha / (alpha + 1)
    in [uS]
gf = 0.00427 [uS]
    in [uS]
alpha = 0.5927
ach_shift = if(ach.ACh <= 0 [mM], 0 [mV],
               -1 [mV] - 9.898 [mV] * (ach.ACh / 1 [mM])^0.618 / ((ach.ACh / 1 [mM])^0.618 + 1.22423e-3))
    in [mV]
iso_shift = if(iso.iso > 0, 7.5 [mV], 0 [mV])
    in [mV]
dot(y) = (inf - y) / tau
    tau = 1 / (0.36 [1/mV/s] * (V + 148.8 [mV] - ach_shift - iso_shift) / (exp(0.066 [1/mV] * (V + 148.8 [mV] - ach_shift - iso_shift)) - 1) + 0.1 [1/mV/s] * (V + 87.3 [mV] - ach_shift - iso_shift) / (1 - exp(-0.2 [1/mV] * (V + 87.3 [mV] - ach_shift - iso_shift)))) - 0.054 [s]
        in [s]
    inf = if(V < -(80 [mV] - ach_shift - iso_shift),
             0.01329 + 0.99921 / (1 + exp((V + 97.134 [mV] - ach_shift - iso_shift) / 8.1752 [mV])),
             0.0002501 * exp((V - ach_shift - iso_shift) / -12.861 [mV]))

#
# Sodium-potassium pump current
#
[inak]
use membrane.V
use extra.K_o, sodium.Na_i, rev.E_Na
Iso_increase = piecewise(iso.iso > 0, 1.2, 1)
Km_Kp = 1.4 [mM]
    in [mM]
Km_Nap = 14 [mM]
    in [mM]
INaK = Iso_increase * INaK_max * (1 + (Km_Kp / K_o)^1.2)^(-1) * (1 + (Km_Nap / Na_i)^1.3)^(-1) * (1 + exp(-(V - E_Na + 110 [mV]) / 20 [mV]))^(-1)
    in [nA]
INaK_max = 0.08105 [nA]
    in [nA]

#
# Sodium-calcium exchanger
#
[inaca]
use membrane.V, phys.RTF
use extra.Na_o, sodium.Na_i
use extra.Ca_o, calcium.Ca_sub
K_NaCa = 3.343 [nA]
    in [nA]
INaCa = K_NaCa * r
    in [nA]
# Cycle rate
r = (x2 * k21 - x1 * k12) / (x1 + x2 + x3 + x4)
x1 = k41 * k34 * (k23 + k21) + k21 * k32 * (k43 + k41)
x2 = k32 * k43 * (k14 + k12) + k41 * k12 * (k34 + k32)
x3 = k14 * k43 * (k23 + k21) + k12 * k23 * (k43 + k41)
x4 = k23 * k34 * (k14 + k12) + k14 * k21 * (k34 + k32)
# Rates
di = 1 + Ca_sub / Kci * (1 + exp(-Qci * V / phys.RTF) + Na_i / Kcni) + Na_i / K1ni * (1 + Na_i / K2ni * (1 + Na_i / K3ni))
do = 1 + Ca_o / Kco * (1 + exp(Qco * V / phys.RTF)) + Na_o / K1no * (1 + Na_o / K2no * (1 + Na_o / K3no))
k12 = Ca_sub / Kci * exp(-Qci * V / RTF) / di
k21 = Ca_o / Kco * exp(Qco * V / RTF) / do
k23 = Na_o / K1no * Na_o / K2no * (1 + Na_o / K3no) * exp(-Qn * V / (2 * RTF)) / do
k32 = exp(Qn * V / (2 * phys.RTF))
k34 = Na_o / (K3no + Na_o)
k43 = Na_i / (K3ni + Na_i)
k14 = Na_i / K1ni * Na_i / K2ni * (1 + Na_i / K3ni) * exp(Qn * V / (2 * RTF)) / di
k41 = exp(-Qn * V / (2 * RTF))
Qci = 0.1369
Qco = 0
Qn = 0.4315
K1ni = 395.3 [mM]
    in [mM]
K1no = 1628 [mM]
    in [mM]
K2ni = 2.289 [mM]
    in [mM]
K2no = 561.4 [mM]
    in [mM]
K3ni = 26.44 [mM]
    in [mM]
K3no = 4.663 [mM]
    in [mM]
Kci = 0.0207 [mM]
    in [mM]
Kcni = 26.44 [mM]
    in [mM]
Kco = 3.663 [mM]
    in [mM]

#
# Calcium release from the SR
#
[carel]
use calcium.Ca_sub, calcium.Ca_jsr
JSRCarel = ks * o * (Ca_jsr - Ca_sub)
    in [mM/s]
ks = 1.480410851e8 [1/s]
    in [1/s]
kom = 660 [1/s]
    in [1/s]
kim = 5 [1/s]
    in [1/s]
kiCa = 500 [1/mM/s]
    in [1/mM/s]
koCa = 10000 [1/mM^2/s]
    in [1/mM^2/s]
ec50_SR = 0.45 [mM]
    in [mM]
MaxSR = 15
MinSR = 1
kCaSR = MaxSR - (MaxSR - MinSR) / (1 + (ec50_SR / Ca_jsr)^2.5)
koSRCa = koCa / kCaSR
    in [1/mM^2/s]
kiSRCa = kiCa * kCaSR
    in [1/mM/s]
dot(r) = kim * ri - kiSRCa * Ca_sub * r - (koSRCa * Ca_sub^2 * r - kom * o)
dot(o) = koSRCa * Ca_sub^2 * r - kom * o - (kiSRCa * Ca_sub * o - kim * i)
dot(i) = kiSRCa * Ca_sub * o - kim * i - (kom * i - koSRCa * Ca_sub^2 * ri)
dot(ri) = kom * i - koSRCa * Ca_sub^2 * ri - (kim * ri - kiSRCa * Ca_sub * r)

#
# Calcium buffering
#
[cabuf]
CM_tot = 0.045 [mM]
    in [mM]
CQ_tot = 10 [mM]
    in [mM]
Mgi = 2.5 [mM]
    in [mM]
TC_tot = 0.031 [mM]
    in [mM]
TMC_tot = 0.062 [mM]
    in [mM]
delta_fCMi = kf_CM * calcium.Ca_i * (1 - fCMi) - kb_CM * fCMi
    in [1/s]
delta_fCMs = kf_CM * calcium.Ca_sub * (1 - fCMs) - kb_CM * fCMs
    in [1/s]
delta_fCQ = kf_CQ * calcium.Ca_jsr * (1 - fCQ) - kb_CQ * fCQ
    in [1/s]
delta_fTC = kf_TC * calcium.Ca_i * (1 - fTC) - kb_TC * fTC
    in [1/s]
delta_fTMC = kf_TMC * calcium.Ca_i * (1 - (fTMC + fTMM)) - kb_TMC * fTMC
    in [1/s]
delta_fTMM = kf_TMM * Mgi * (1 - (fTMC + fTMM)) - kb_TMM * fTMM
    in [1/s]
dot(fCMi) = delta_fCMi
dot(fCMs) = delta_fCMs
dot(fCQ) = delta_fCQ
dot(fTC) = delta_fTC
dot(fTMC) = delta_fTMC
dot(fTMM) = delta_fTMM
kb_CM = 542 [1/s]
    in [1/s]
kb_CQ = 445 [1/s]
    in [1/s]
kb_TC = 446 [1/s]
    in [1/s]
kb_TMC = 7.51 [1/s]
    in [1/s]
kb_TMM = 751 [1/s]
    in [1/s]
kf_CM = 1642000 [1/mM/s]
    in [1/mM/s]
kf_CQ = 175.4 [1/mM/s]
    in [1/mM/s]
kf_TC = 88800 [1/mM/s]
    in [1/mM/s]
kf_TMC = 227700 [1/mM/s]
    in [1/mM/s]
kf_TMM = 2277 [1/mM/s]
    in [1/mM/s]

#
# Intracellular calcium fluxes
#
[caflux]
JCa_dif = (calcium.Ca_sub - calcium.Ca_i) / 5.469e-5 [s]
    in [mM/s]
Jtr = (calcium.Ca_nsr - calcium.Ca_jsr) / 0.04 [s]
    in [mM/s]
Jup = P_up / (1 + exp((-calcium.Ca_i + K_up) / 5e-5 [mM]))
    in [mM/s]
K_up =  2.86113e-4 [mM]
    in [mM]
P_up = P_up_basal * (1 - b_up)
    in [mM/s]
P_up_basal = 5 [mM/s]
    in [mM/s]
b_up = piecewise(iso.iso > 0, -0.25, ach.ACh > 0 [mM], 0.7 * ach.ACh / (9e-5 [mM] + ach.ACh), 0)

#
# Intracellular calcium concentrations
#
[calcium]
use cell.vsub, cell.vi, cell.vjsr, cell.vnsr
dot(Ca_jsr) = caflux.Jtr - (carel.JSRCarel + cabuf.CQ_tot * cabuf.delta_fCQ)
    in [mM]
dot(Ca_nsr) = caflux.Jup - caflux.Jtr * vjsr / vnsr
    in [mM]
dot(Ca_sub) = carel.JSRCarel * vjsr / vsub - ((ical.IsiCa + icat.ICaT - 2 * inaca.INaCa) / (2 * phys.F * vsub) + caflux.JCa_dif + cabuf.CM_tot * cabuf.delta_fCMs)
    in [mM]
dot(Ca_i) = (caflux.JCa_dif * vsub - caflux.Jup * vnsr) / vi - (cabuf.CM_tot * cabuf.delta_fCMi + cabuf.TC_tot * cabuf.delta_fTC + cabuf.TMC_tot * cabuf.delta_fTMC)
    in [mM]

#
# Intracellular sodium concentration
#
[sodium]
dot(Na_i) = (1 - clamp) * -INa_tot / ((cell.vi + cell.vsub) * phys.F)
    in [mM]
INa_tot = ina.INa + if.IfNa + ical.IsiNa + 3 * inak.INaK + 3 * inaca.INaCa
    in [nA]
clamp = 1

#
# Intracellular potassium concentration
#
[potassium]
K_i = 140 [mM]
    in [mM]

[[protocol]]
# This model is self-excitatory

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

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

# Run simulation
d = s.run(3)

# Display the results
plt.figure()
plt.plot(d.time(), d['membrane.V'])
plt.show()