ni-2017.mmt

[[model]]
name: ni-2017
version: 20240904
mmt_authors: Michael Clerx
display_name: Ni et al., 2017
desc: """
    The 2017 model of the human atrial AP by Ni et al. [1].

    This implementation is based on the code for the model by [1], which is an
    update of a model by Colman et al. (2013, 2017) [2, 3], which in turn was
    based on the 1998 model by Courtemanche et al. [4]. The model includes a
    simplified version of the calcium subspaces used by Koivumaki et al. [5].

    The state derivatives from this implementation were checked against output
    from the original code.

    The model contained a mixture of volume and time units (pL, inherited from
    [5] and um^3, inherited from [4]). This was adjusted so that all rates are
    in 1/ms and all volumes in um^3. In doing so, a scaling mistake in the
    equations for the SR calcium concentrations (resulting in the dynamics
    being sped up). This was left uncorrected, so that this implementation
    matches the published code.

    In this Myokit implementation, the stimulus was set to 0.5 ms and
    approximately twice the threshold value for depolarisation.

    [1] Ni, H., Whittaker, D. G., Wang, W., Giles, W. R., Narayan, S. M., &
        Zhang, H. (2017). Synergistic Anti-arrhythmic Effects in Human Atria
        with Combined Use of Sodium Blockers and Acacetin. Frontiers in
        Physiology, 8, 946.
        https://doi.org/10.3389/fphys.2017.00946

    [2] Colman, M. A., Aslanidi, O. V., Kharche, S., Boyett, M. R., Garratt,
        C., Hancox, J. C., & Zhang, H. (2013). Pro-arrhythmogenic effects of
        atrial fibrillation-induced electrical remodelling: insights from the
        three-dimensional virtual human atria. Journal of Physiology, 591(17),
        4249-4272.
        https://doi.org/10.1113/jphysiol.2013.254987

    [3] Colman, M. A., Ni, H., Liang, B., Schmitt, N., & Zhang, H. (2017). In
        silico assessment of genetic variation in KCNA5 reveals multiple
        mechanisms of human atrial arrhythmogenesis. PLOS Computational
        Biology, 13(6), e1005587.
        https://doi.org/10.1371/journal.pcbi.1005587

    [4] Courtemanche, M., Ramirez, R. J., & Nattel, S. (1998). Ionic mechanisms
        underlying human atrial action potential properties: insights from a
        mathematical model. American Journal of Physiology. Heart and
        Circulatory Physiology, 275(1), H301-H321.
        https://doi.org/10.1152/ajpheart.1998.275.1.H301

    [5] Koivumaki, J. T, Korhonen, T., & Tavi P. (2011). Impact of Sarcoplasmic
        Reticulum Calcium Release on Calcium Dynamics and Action Potential
        Morphology in Human Atrial Myocytes: A Computational Study. PLOS
        Computational Biology, 7(1), e1001067.
        https://doi.org/10.1371/journal.pcbi.1001067

    Original comments follow below:

    ---------------------------------------------------------------------------

    Implementation of IKur and INa state-dependent drug blockade. (using
    guarded-receptor model)

    This is the updatd CNZ model, update from the Colman et al. model 2013
    Journal of Physiology Paper,

    Updates:  New IKur Model,
    Modifications to the conductances of several key ionic currents.

    Modifications to the Calcium dynamics model. change in the Refractory
    Parameters to make it faster, and thus remove the alternans seen at 400 ms
    of BCL.

    RMP: around -76 mV.
"""
# Initial values
membrane.V      = -77.13255836
ina.m           = 0.005631819916
ina.h           = 0.9168420281
ina.j           = 0.9380183564
ical.d          = 0.0002267161277
ical.f          = 0.9354212881
ikr.xr          = 0.00157418133
iks.xs          = 0.02225979641
sodium.Na_i     = 10.30397012
calcium.Ca_i    = 0.0001403133065
potassium.K_i   = 131.867138
ical.fca        = 0.7270823666
ito.r           = 0.01223706011
ito.s           = 0.8849139842
calcium.Ca_ss   = 0.0001313595105
calcium.Ca_SRi  = 0.9892411621
calcium.Ca_SRss = 0.9779168037
serca.SERCACa   = 0.00958584702
serca.SERCACass = 0.009386941118
ryr.o_ss        = 0.000456694441
ryr.c_ss        = 0.9571976502
ryr.a_ss        = 0.1419271573
ryr.o_i         = 0.0003667999273
ryr.c_i         = 0.9790238698
ryr.a_i         = 0.2977219443
ikur.a          = 0.0002417881801
ikur.i          = 0.9517278864
ina.BA          = 0
ina.BI          = 0
ikur.BO         = 2.355477517e-16
ikur.BC         = 3.062722822e-15


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

#
# Membrane potential
#
[membrane]
use stimulus.i_stim
i_ion = (
        + ina.INa
        + ito.Ito
        + ikur.IKur
        + ikr.IKr
        + iks.IKs
        + ical.ICaL
        + ik1.IK1
        + ibna.IbNa
        + ibca.IbCa
        + inaca.INaCa
        + inak.INaK
        + icap.ICap
    )
    in [pA]
    label cellular_current
dot(V) = -(i_ion + i_stim + i_diff) / cell.Cm
    in [mV]
    label membrane_potential
i_diff = 0 [pA]
    in [pA]
    bind diffusion_current

#
# Stimulus current
#
[stimulus]
i_stim = cell.Cm * engine.pace * amplitude
    in [pA]
amplitude = -36 [A/F] * 2
    in [pA/pF]

#
# AF mode
#
[af]
af = 0
    desc: Set to 1 to enable AF mode

#
# Physical constants
#
[phys]
T = 310 [K] # 37 Celcius
    in [K]
F = 96.4867 [C/mmol]
    in [C/mmol]
R = 8.3143 [J/mol/K]
    in [J/mol/K]
FRT = F / R / T
    in [1/mV]

#
# Cell geometry
#
[cell]
Cm = 100 [pF]
    in [pF]
    label membrane_capacitance

# Diffusion compartment volumes

# Bulk cytosol volume used by Na (and potentially K)
# This value is from Courtemanche et al.
Vi = 13668 [um^3]
    in [um^3]

# Junctional subspace / subsarcolemmal subspace.
# This is from Koivumaki et al., and represents a space just below the membrane
# The value is multiplied by 2, as this cell is 100pF (vs 50pF in Koivumaki)
Vss = 2 * 49.9232 [um^3]
    in [um^3]

# Volume of the non-junctional space that Cai can diffuse in
Vnj = 6 * 2531 [um^3]          # Koviumaki value times 6
    in [um^3]                  # The name of this variable comes from Koivumaki

# Volume of the SR interacting with the bulk and junctional spaces
VSRi = 2 * 57 [um^3]           # Koivumaki value times 2
    in [um^3]                  # VSRi and VSRss are also koivumaki variables
VSRss = 2 * 80 [um^3]          # Koivumaki value times 2
    in [um^3]

# Ca diffusion grid (from Koivumaki et al.)
dx = 1.625 [um]
    in [um]
dx2 = dx * dx
    in [um^2]
# Surface area between the junctional (SS) and non junctional (bulk) spaces
Aj_nj = 2492.324412 [um^2]  # Koivumaki value
    in [um^2]
# Diffusion distance from center to center of junct to first njunct
xj_nj = 0.822500 [um]       # Koivumaki value
    in [um]

#
# External ion concentrations
#
[ion]
Na_o = 140 [mM]
    in [mM]
Ca_o = 1.8 [mM]
    in [mM]
K_o = 5.4 [mM]
    in [mM]

#
# Reversal potentials
#
[nernst]
ENa = 26.71 [mV] * log(ion.Na_o / sodium.Na_i)
    in [mV]
EK = 26.71 [mV] * log(ion.K_o / potassium.K_i)
    in [mV]
ECa = 13.35 [mV] * log(ion.Ca_o / calcium.Ca_i)
    in [mV]

#
# Fast sodium current
#
[ina]
use membrane.V, cell.Cm

# Drug binding
drug_concen = 0 [uM]
    in [uM]
drug_Ka = 0.1 [1/ms/uM] / 1000
    in [1/ms/uM]
drug_Ki = 0.1 [1/ms/uM] / 1000
    in [1/ms/uM]
drug_La = 0.1 [1/ms]
    in [1/ms]
drug_Li = 0.01 [1/ms]
    in [1/ms]

dot(BA) = drug_Ka * drug_concen * m ^3 * h * j * (1 - BA - BI) - drug_La * BA
    desc: blockade to Open gate
dot(BI) = drug_Ki * drug_concen * (1 - h) * (1 - BA - BI) - drug_Li * BI
    desc: blockade to closed states

# Gating
dot(m) = alpha * (1 - m) - beta * m
    alpha = if(abs(V - shift + 47.13 [mV]) < 1e-10 [mV],
               3.2 [1/ms],
               0.32 [1/ms/mV] * (V - shift + 47.13 [mV]) / (1 - exp(-0.1 [1/mV] * (V - shift + 47.13 [mV]))))
        in [1/ms]
    beta = 0.08 [1/ms] * exp(-(V - shift) / 11 [mV])
        in [1/ms]
    shift = if(af.af == 1, 1.6 [mV], 0 [mV])
        in [mV]
dot(h) = alpha * (1 - h) - beta * h
    alpha = if(V >= -40 [mV],
               0 [1/ms],
               0.135 [1/ms] * exp((V + 80 [mV]) / -6.8 [mV]))
        in [1/ms]
    beta = if(V >= -40 [mV],
              1 / (0.13 [ms] * (1 + exp((V + 10.66 [mV]) / -11.1 [mV]))),
              3.56 [1/ms] * exp(0.079 [1/mV] * V) + 3.1e5 [1/ms] * exp(0.35 [1/mV] * V))
        in [1/ms]
dot(j) = alpha * (1 - j) - beta * j
    alpha = if(V >= -40 [mV],
               0 [1/ms],
               (-1.2714e5 [1/ms/mV] * exp(0.2444 [1/mV] * V) - 3.474e-5 [1/ms/mV] * exp(-0.04391 [1/mV] * V)) * (V + 37.78 [mV]) / (1 + exp(0.311 [1/mV] * (V + 79.23 [mV]))))
        in [1/ms]
    beta = if(V >= -40 [mV],
              0.3 [1/ms] * exp(-2.535e-7 [1/mV] * V) / (1 + exp(-0.1 [1/mV] * (V + 32 [mV]))),
              0.1212 [1/ms] * exp(-0.01052 [1/mV] * V) / (1 + exp(-0.1378 [1/mV] * (V + 40.14 [mV]))))
        in [1/ms]
gNa = 7.8 [nS/pF]
    in [nS/pF]
INa = gNa * m ^3 * h * j * (V - nernst.ENa) * Cm * (1 - BA - BI)
    in [pA]

#
# Rapidly activating K current
#
[ikr]
use membrane.V, cell.Cm
use acacetin.fKr
dot(xr) = (inf - xr) / tau
    a = if(abs(V + 14.1 [mV]) < 1e-10 [mV],
           0.0015 [1/ms],
           0.0003 [1/ms/mV] * (V + 14.1 [mV]) / (1 - exp((V + 14.1 [mV]) / -5 [mV])))
        in [1/ms]
    b = if(abs(V - 3.3328 [mV]) < 1e-10 [mV],
           3.7836118e-4 [1/ms],
           0.000073898 [1/ms/mV] * (V - 3.3328 [mV]) / (exp((V - 3.3328 [mV]) / 5.1237 [mV]) - 1))
        in [1/ms]
    tau = 1 / (a + b)
        in [ms]
    inf = 1 / (1 + exp((V + 14.1 [mV]) / -6.5 [mV]))
gKr = 0.8 * 0.029411765 [nS/pF]
    in [nS/pF]
    label g_Kr
IKr = Cm * fKr * gKr * xr * (V - nernst.EK) / (1 + exp((V + 15 [mV]) / 22.4 [mV]))
    in [pA]

#
# Slowly activating K current
#
[iks]
use membrane.V, cell.Cm
use acacetin.fKs
dot(xs) = (inf - xs) / tau
    a = if(abs(V - 19.9 [mV]) < 1e-10 [mV],
           0.00068 [1/ms],
           0.00004 [1/ms/mV] * (V - 19.9 [mV]) / (1 - exp((V - 19.9 [mV]) / -17 [mV])))
        in [1/ms]
    b = if(abs(V - 19.9 [mV]) < 1e-10 [mV],
           0.000315 [1/ms],
           0.000035 [1/ms/mV] * (V - 19.9 [mV]) / (exp((V - 19.9 [mV]) / 9 [mV]) - 1))
        in [1/ms]
    tau = 0.5 / (a + b)
        in [ms]
    inf = sqrt(1 / (1 + exp((V - 19.9 [mV]) / -12.7 [mV])))
gKs = 0.8 * 0.12941176 [nS/pF] * if(af.af == 1, 2, 1)
    in [nS/pF]
IKs = Cm * fKs * gKs * xs * xs * (V - nernst.EK)
    in [pA]

#
# L-type calcium current
#
[ical]
use membrane.V, cell.Cm
dot(d) = (inf - d) / tau
    inf = 1 / (1 + exp((V + 10 [mV]) / -8 [mV]))
    a = 1 [ms] / (1 + exp((V + 10 [mV]) / -6.24 [mV]))
        in [ms]
    tau = if(abs(V + 10 [mV]) < 1e-10 [mV],
             a * 4.579,
             a * (1 - exp((V + 10 [mV]) / -6.24 [mV])) / (0.035 [1/mV] * (V + 10 [mV])))
        in [ms]
dot(f) = (inf - f) / tau
    inf = exp(-(V + 28 [mV]) / 6.9 [mV]) / (1 + exp(-(V + 28 [mV]) / 6.9 [mV]))
    tau = 1.5 [ms] * 2 * 3 / (0.0197 * exp(-(0.0337 [1/mV])^2 * (V + 10 [mV])^2) + 0.02)
        in [ms]
dot(fca) = (inf - fca) / tau
    inf = 1 / (1 + (calcium.Ca_ss / 0.00035 [mM]))
    tau = 2 [ms]
        in [ms]
gCaL = 0.1294 [nS/pF] * 0.75 * if(af.af == 1, 0.35, 1)
    in [nS/pF]
ErL = 65.0 [mV]
    in [mV]
ICaL = 1.333333 * Cm * gCaL * d * f * fca * (V - ErL)
    in [pA]

#
# Transient outward current
#
[ito]
use membrane.V, cell.Cm
use acacetin.fto
dot(r) = (inf - r) / tau
    inf = 1 / (1 + exp((V - shift - 1 [mV]) / -11 [mV]))
    tau = 3.5 [ms] * exp(-((V - shift) / 30 [mV]) * 2) + 1.5 [ms]
        in [ms]
    shift = if(af.af == 1, 16 [mV], 0 [mV])
        in [mV]
dot(s) = (inf - s) / 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) + 14.14 [ms]
        in [ms]
gto = 0.75471 * 0.1962 [nS/pF] * if(af.af == 1, 0.35, 1)
    in [nS/pF]
Ito = 1.05 * fto * Cm * gto * r * s * (V - nernst.EK)
    in [pA]

#
# Acacetin concentration
#
[acacetin]
conc = 1e-14 [uM]
    in [uM]
Aca_IC50_Ito = 9.3 [uM]
    in [uM]
Aca_IC50_IKr = 32.4 [uM]
    in [uM]
Aca_IC50_IKs = 81.4 [uM]
    in [uM]
Aca_hill_Ito = 0.9
Aca_hill_IKr = 0.9
Aca_hill_IKs = 0.8
fto =  1 - 1 / (1 + (Aca_IC50_Ito / conc)^Aca_hill_Ito)
fKr =  1 - 1 / (1 + (Aca_IC50_IKr / conc)^Aca_hill_IKr)
fKs =  1 - 1 / (1 + (Aca_IC50_IKs / conc)^Aca_hill_IKs)

#
# Ultra rapid delayed rectifier outward K Current
#
[ikur]
use membrane.V, cell.Cm
use phys.FRT

# Drug binding
KO  =  0.00019445959 [1/ms/uM]
    in [1/ms/uM]
LO  =  0.00029099869 [1/ms]
    in [1/ms]
ZKO = -0.25731953173
ZLO = -0.01280158014
KC  =  0.00247871659 [1/ms/uM]
    in [1/ms/uM]
LC  =  0.00028636715 [1/ms]
    in [1/ms]
ZKC = -0.32702357338
ZLC =  0.83707302269

dot(BO) = KO * exp(-ZKO * V * FRT) * acacetin.conc * i * a * (1 - BO - BC) - LO * BO * exp(-ZLO * V * phys.FRT)
dot(BC) = KC * exp(-ZKC * V * FRT) * acacetin.conc * i * (1 - a) * (1 - BO - BC) - LC * BC * exp(-ZLC * V * phys.FRT)

# Gating
K_Q10 = 3.0
dot(a) = (inf - a) / tau
    inf = 1 / (1 + exp(-(V + 5.52 [mV]) / 8.6 [mV]))
    tau = (
            (45.6666746826 / (1 + exp((V + 11.2306497073 [mV]) / 11.5254705962 [mV])) + 4.26753514993) *
            (0.262186042981 / (1 + exp((V + 35.8658312707 [mV]) / -3.87510627762 [mV])) + 0.291755017928)
          ) * 1 [ms] / K_Q10
        in [ms]
dot(i) = (inf - i) / tau
    inf = 0.52424 / (1 + exp((V + 15.1142 [mV]) / 7.567021 [mV])) + 0.4580778
    tau = (2328 [ms] / (1 + exp((V - 9.435 [mV]) / 3.5827 [mV])) + 1739.139 [ms]) / K_Q10
        in [ms]

gKur = 0.006398 [nS/pF] * 0.9 * if(af.af == 1, 0.5, 1)
    in [nS/pF]
    label g_Kur
IKur = Cm * gKur * (4.5128 + 1.899769 / (1 + exp((V - 20.5232 [mV]) / -8.26597 [mV]))) * (1 - BO - BC) * a * i * (V - nernst.EK)
    in [pA]

#
# Time-independent K current
#
[ik1]
use membrane.V, cell.Cm
gK1 = 0.09 [nS/pF] * if(af.af == 1, 2.1, 1)
    in [nS/pF]
IK1 = Cm * gK1 * (V - nernst.EK) / (1 + exp(0.07 [1/mV] * (V + 80 [mV])))
    in [pA]

#
# Background calcium current
#
[ibca]
use membrane.V, cell.Cm
gbCa = 1.4 * 0.001131 [nS/pF]
    in [nS/pF]
IbCa = 1.0 * Cm * gbCa * (V - nernst.ECa)
    in [pA]

#
# Background sodium current
#
[ibna]
use membrane.V, cell.Cm
gbNa = 0.8 * 0.0006744375 [nS/pF]
    in [nS/pF]
IbNa = 1.7 * Cm * gbNa * (V - nernst.ENa)
    in [pA]

#
# Sarcolemmal Ca pump current
#
[icap]
use cell.Cm
use calcium.Ca_ss
icapbar = 0.275 [A/F]
    in [A/F]
kmcap = 0.0005 [mM]
    in [mM]
ICap = 1.26 * Cm * icapbar * Ca_ss / (Ca_ss + kmcap)
    in [pA]

#
# Sodium-calcium exchanger current
#
[inaca]
use membrane.V
use ion.Na_o, ion.Ca_o
knacalr = 1600 [A/F]
    in [A/F]
kmnalr = 87.5 [mM]
    in [mM]
kmcalr = 1.38 [mM]
    in [mM]
ksatlr = 0.1
gammalr = 0.35
scaling = if(af.af == 1, 1.4, 1)
INaCa = (1.4 * scaling * cell.Cm * knacalr / (kmnalr^3 + Na_o^3) / (kmcalr + Ca_o)
         / (1 + ksatlr * exp((gammalr - 1) * V * phys.FRT))
         * (sodium.Na_i^3 * Ca_o * exp(V * gammalr * phys.FRT) - Na_o^3 * calcium.Ca_ss *  exp(V * (gammalr - 1) * phys.FRT))
        )
    in [pA]

#
# Sodium-potasisum pump
#
[inak]
use membrane.V, cell.Cm, phys.FRT
use ion.Na_o, ion.K_o, sodium.Na_i
sigma = (exp(Na_o / 67.3 [mM]) - 1) / 7
fnak = 1 / (1 + 0.1245 * exp(-0.1 * V * FRT) + 0.0365 * sigma * exp(-V * FRT))
kmko = 1.5 [mM]
    in [mM]
kmnai = 10 [mM]
    in [mM]
I_bar = 1.4 * 0.59933874 [A/F]
    in [A/F]
INaK = 1.28 * Cm * I_bar * fnak * K_o / (K_o + kmko) / (1 + (kmnai / Na_i)^4)
    in [pA]

#
# Intercellular sodium concentration
#
[sodium]
dot(Na_i) = (-3 * inak.INaK - 3 * inaca.INaCa - ibna.IbNa - ina.INa) / (phys.F * cell.Vi)
    in [mM]

#
# Intercellular potassium concentration
#
[potassium]
#dot(K_i) = (2 * inak.INaK - ik1.IK1 - ito.Ito - ikur.IKur - ikr.IKr - iks.IKs) / (phys.F * cell.Vi)
# Original comment: K+ is drifting
dot(K_i) = 0 [mM/ms]
    in [mM]

#
# Ryanodine release currents
# This part is adapted from Koivumaki et al.
#
[ryr]
use calcium.Ca_i, calcium.Ca_ss, calcium.Ca_SRi, calcium.Ca_SRss
fRyr = if(af.af == 1, 2.5, 0)
tau_scaling = if(af.af == 1, 2.7, 1)
dot(o_i) = (inf - o_i) / tau
    inf = 1 - 1 / (1 +  exp((2000 [1/mM] * Ca_i * (1 + fRyr) - (a_i + 0.22)) / 0.03))
    tau = 5 [ms]
        in [ms]
dot(c_i) = (inf - c_i) / tau
    inf = 1 / (1 +  exp((2000 [1/mM] * Ca_i * (1 + fRyr) - (a_i + 0.02)) / 0.01))
    tau = 2 * 15 [ms] * tau_scaling
        in [ms]
dot(a_i) = (inf - a_i) / tau
    inf = 0.505 - 0.427 / (1 + exp((2000 [1/mM] * Ca_i * (1 + fRyr) - 0.29) / 0.082))
    tau = 250 [ms] * tau_scaling
        in [ms]
dot(o_ss) = (inf - o_ss) / tau
    inf = 1 - 1 / (1 +  exp((1000 [1/mM] * Ca_ss * (1 + fRyr) - (a_ss + 0.22)) / 0.03))
    tau = 5 [ms]
        in [ms]
dot(c_ss) = (inf - c_ss) / tau
    inf = 1 / (1 + exp((1000 [1/mM] * Ca_ss * (1 + fRyr) - (a_ss + 0.02)) / 0.01))
    tau = 15 [ms] * tau_scaling
        in [ms]
dot(a_ss) = (inf - a_ss) / tau
    inf = 0.505 - 0.427 / (1 + exp((1000 [1/mM] * Ca_ss * (1 + fRyr) - 0.29) / 0.082))
    tau = 250 [ms] * tau_scaling
        in [ms]

scaling = if(af.af == 1, 3, 1)

# Release from SR into the bulk cytosol
nui = 1e-3 [1/ms] * cell.Vnj
    in [um^3/ms]
RyRSRCai = 1 - 1 / (1 +  exp((Ca_SRi - 0.3 [mM]) / 0.1 [mM]))
Jreli = scaling * nui * o_i * c_i * RyRSRCai * (Ca_SRi - Ca_i)
    in [um^3*mM/ms]

# Release from SR into the junctional subspace
nuss = 0.625 [1/ms] * cell.Vss
    in [um^3/ms]
RyRSRCass = 1 - 1 / (1 + exp((Ca_SRss - 0.3 [mM]) / 0.1 [mM]))
Jrelss = scaling * nuss * o_ss * c_ss * RyRSRCass * (Ca_SRss - Ca_ss)
     in [um^3*mM/ms]

#
# SERCA
# Adapted from Koivumaki et al.
# The fluxes J are fluxes in 1um^3 of volume.
#
[serca]
use cell.Vss, cell.Vnj
k1 = 7500 [1/ms/mM^2]
    desc: 1000^2 * k4
    in [1/ms/mM^2]
k2 = 4.6875e-4 [1/ms]
    desc: k1 * SERCAKmf^2
    in [1/ms]
k3 = 2.314815e-3 [1/ms/mM^2]
    desc: k4 / SERCAKmr^2
    in [1/ms/mM^2]
k4 = 7.5e-3 [1/ms]
    desc: pump rate
    in [1/ms]

scaling = if(af.af == 1, 0.6, 1)
cpumps = 40e-3 [mM]
    in [mM]
J_SERCASR = scaling * 0.75 * (-k3 * calcium.Ca_SRi^2 * (cpumps - SERCACa) + k4 * SERCACa) * Vnj * 2
    in [um^3*mM/ms]
J_SERCASRss = scaling * 0.75 * (-k3 * calcium.Ca_SRss^2 * (cpumps - SERCACass) + k4 * SERCACass) * Vss * 2
    in [um^3*mM/ms]
J_bulkSERCAss = scaling * 0.75 * (k1 * calcium.Ca_ss^2 * (cpumps - SERCACass) - k2 * SERCACass) * Vss * 2
    in [um^3*mM/ms]
J_bulkSERCA = scaling * 0.75 * (k1 * calcium.Ca_i^2 * (cpumps - SERCACa) - k2 * SERCACa) * Vnj * 2
    in [um^3*mM/ms]

# Calcium buffered by SERCA:
dot(SERCACa) = (-J_SERCASR + J_bulkSERCA) / Vnj
    in [mM]
dot(SERCACass) = (-J_SERCASRss + J_bulkSERCAss) / Vss
    in [mM]

#
# Internal calcium concentrations
# Mostly based on Koivumaki
#
[calcium]
use cell.Vss, cell.Vnj, cell.dx2

# Diffusion from junctional to non-junctional
DCa = 0.78 [um^2/ms]
    in [um^2/ms]
Jj_nj = 2.5 * DCa * cell.Aj_nj / cell.xj_nj * (Ca_ss - Ca_i)
    in [um^3*mM/ms]

scaling = if(af.af == 1, 1.5, 1)
kSRleak = 6e-6 [1/ms]
    in [1/ms]
JSRCaleaki = 0.5 * scaling * kSRleak * (Ca_SRi - Ca_i) * Vnj
    in [um^3*mM/ms]
JSRCaleakss = 0.5 * scaling * kSRleak * (Ca_SRss - Ca_ss) * Vss
    in [um^3*mM/ms]

# Ca fluxes in 1 um^3 volume
JCa = -serca.J_bulkSERCA + JSRCaleaki + ryr.Jreli + Jj_nj
    in [um^3*mM/ms]
JCass = -Jj_nj + JSRCaleakss - serca.J_bulkSERCAss + ryr.Jrelss
    in [um^3*mM/ms]
JSRCai = serca.J_SERCASR - JSRCaleaki - ryr.Jreli
    in [um^3*mM/ms]
JSRCaSS = serca.J_SERCASRss - JSRCaleakss - ryr.Jrelss
    in [um^3*mM/ms]

# Calcium concentration in the bulk (non-junctional) space
KdBCa = 2.38e-3 [mM]
    in [mM]
BCa = 24e-3 [mM]
    in [mM]
dot(Ca_i) = JCa / Vnj * beta
    beta = 1 / (1 + BCa * KdBCa / (Ca_i + KdBCa)^2)
    in [mM]

# Calcium stored in the SR
DCaSR = 0.044 [um^2/ms]
    in [um^2/ms]
CSQN = 6.7 [mM]
    in [mM]
KdCSQN = 0.8 [mM]
    in [mM]
dot(Ca_SRi) = 1000 * beta * (DCaSR * ((Ca_SRss - 2 * Ca_SRi + Ca_SRi) / dx2 + (Ca_SRi - Ca_SRss) / (2 * 3 * dx2)) + JSRCai / cell.VSRi)
    beta = 1 / (1 + CSQN * KdCSQN / (Ca_SRi + KdCSQN)^2)
    in [mM]
dot(Ca_SRss) = 1000 * beta * (DCaSR * ((Ca_SRss - 2 * Ca_SRss + Ca_SRi) / dx2 + (Ca_SRss - Ca_SRi) / (2 * 4 * dx2)) + JSRCaSS / cell.VSRss)
    beta = 1 / (1 + CSQN * KdCSQN / (Ca_SRss + KdCSQN)^2)
    in [mM]
# N.B.: The factor 1000 in the above equations was a mistake, but is left in
# here to match the original code.

# SS Ca
SLlow = 165 [mM]
    in [mM]
SLhigh = 13 [mM]
    in [mM]
KdSLlow = 1.1 [mM]
    in [mM]
KdSLhigh = 13e-3 [mM]
    in [mM]
I_tot = -ical.ICaL - ibca.IbCa - icap.ICap + 2 * inaca.INaCa
    in [pA]
dot(Ca_ss) = beta * (JCass / Vss + I_tot / (2 * Vss * phys.F))
    beta = 1 / (1 + SLlow * KdSLlow / (Ca_ss + KdSLlow)^2 + SLhigh * KdSLhigh / (Ca_ss + KdSLhigh)^2 + BCa * KdBCa / (Ca_ss + KdBCa)^2)
    in [mM]

[[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(100 * 1000)
d = s.run(1000)

# Display the results
plt.figure()
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential (mV)')
plt.plot(d.time(), d['membrane.V'])
plt.show()