carro-2011.mmt

[[model]]
name: carro-2011
version: 20240904
mmt_authors: Michael Clerx
display_name: Carro et al., 2011
desc: """
    The 2011 model of the human ventricular AP by Carro et al. [1], based
    on the model by Grandi et al. [2].

    This implementation is based on the CellML code provided by the study
    authors [3]. 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] Carro, J., Rodríguez, J. F., Laguna, P., & Pueyo, E. (2011). A human
        ventricular cell model for investigation of cardiac arrhythmias under
        hyperkalaemic conditions. Philosophical Transactions of the Royal
        Society A: Mathematical, Physical and Engineering Sciences, 369(1954),
        4205-4232.
        https://doi.org/10.1098/rsta.2011.0127

    [2] Grandi, E., Pasqualini, F. S., & Bers, D. M. (2010). A novel
        computational model of the human ventricular action potential and Ca
        transient. Journal of Molecular and Cellular Cardiology, 48, 112-121.
        https://doi.org/10.1016/j.yjmcc.2009.09.019

    [3] https://models.physiomeproject.org/workspace/66b
        Accessed on 2024-09-03

"""
# Initial values
membrane.V    = -84.13368
ina.m         = 2.163678e-3
ina.h         = 0.713497
ina.j         = 0.7128671
ical.d        = 1.871177e-6
ical.f        = 0.9804391
ical.f2       = 0.999401
ical.fCa_jn   = 2.847118e-2
ical.fCa_sl   = 1.692189e-2
ito.xs        = 3.584727e-4
ito.ys        = 0.8087629
ito.xf        = 3.584625e-4
ito.yf        = 0.9999976
ikr.x         = 1.516232e-2
iks.x         = 3.549354e-3
ryr.r         = 0.8886338
ryr.o         = 1.126209e-6
ryr.i         = 1.411382e-7
nabss.NaB_jn  = 3.796195
nabss.NaB_sl  = 0.8283308
cabct.TnCL    = 9.757237e-3
cabct.TnCHCa  = 0.1225914
cabct.TnCHMg  = 8.12201e-3
cabct.CaM     = 3.267494e-4
cabct.MyoCa   = 2.520383e-3
cabct.MyoMg   = 0.1369529
cabct.SRB     = 2.373753e-3
cabss.SLL_jn  = 8.563314e-3
cabss.SLL_sl  = 1.097424e-2
cabss.SLH_jn  = 8.053908e-2
cabss.SLH_sl  = 0.1235381
cabsr.Csqn    = 1.258048
cabsr.Ca_sr   = 0.6093596
sodium.Na_jn  = 10.07931
sodium.Na_sl  = 10.0781
sodium.Na_i   = 10.07825
calcium.Ca_jn = 2.038197e-4
calcium.Ca_sl = 1.184305e-4
calcium.Ca_i  = 9.658067e-5

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

#
# Cell type
#
[mode]
epi = 1
    desc: 1 For epicardial cell, 0 for endocardial cell

#
# Membrane potential
#
[membrane]
use stimulus.i_stim
dot(V) = -(i_ion + i_stim)
    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 + potassium.I_K_tot
    in [A/F]
    label cellular_current

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

#
# Physical constants
#
# Unchanged from [2].
#
[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]

#
# Cell geometry
#
# Unchanged from [2].
#
[cell]
C = 1.381e-10 [F]
    in [F]
L = 100 [um]
    in [um]
R = 10.25 [um]
    in [um]
pi = 3.14159265358979312
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.000539 * 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]
K_i = 138 [mM]  # Increased from 120 in [2]
    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
#
# Unchanged from [2].
#
[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
#
# Rescaled from [2], but otherwise unchanged.
#
[ina]
use membrane.V
INa = INa_jn + INa_sl
    in [A/F]
INa_jn = cell.fjn * gNa * m^3 * h * j * (V - rev.ENa_jn)
    in [A/F]
INa_sl = cell.fsl * gNa * m^3 * h * j * (V - rev.ENa_sl)
    in [A/F]
gNa = 18.86 [mS/uF]
    in [mS/uF]
dot(m) = (inf - m) / tau
    inf = 1 / (1 + exp(-(56.86 [mV] + V) / 9.03 [mV]))^2
    tau = 0.1292 [ms] * exp(-((V + 45.79 [mV]) / 15.54 [mV])^2) + 0.06487 [ms] * exp(-((V - 4.823 [mV]) / 51.12 [mV])^2)
        in [ms]
dot(h) = (inf - h) / tau
    a = if(V >= -40 [mV], 0 [1/ms],
            0.057 [1/ms] * exp(-(V + 80 [mV]) / 6.8 [mV]))
        in [1/ms]
    b = if(V >= -40 [mV],
            5.9231 [1/ms] / (1 + exp(-(V + 10.66 [mV]) / 11.1 [mV])),
            2.7 [1/ms] * exp(0.079 [1/mV] * V) + 3.1e5 [1/ms] * exp(0.3485 [1/mV] * V))
        in [1/ms]
    tau = 1 / (a + b)
        in [ms]
    inf = 1 / (1 + exp((V + 71.55 [mV]) / 7.43 [mV]))^2
dot(j) = (inf - j) / tau
    a = if(V >= -40 [mV], 0 [1/ms], (
                (-2.5428e4 [1/mV/ms] * exp(0.2444 [1/mV] * V) - 6.948e-6 [1/mV/ms] * exp(-0.04391 [1/mV] * V))
                * (V + 37.78 [mV]) / (1 + exp(0.311 [1/mV] * (V + 79.23 [mV])))
            ))
        in [1/ms]
    b = if(V >= -40 [mV],
            0.6 [1/ms] * exp(0.057 [1/mV] * V) / (1 + exp(-0.1 [1/mV] * (V + 32 [mV]))),
            0.02424 [1/ms] * exp(-0.01052 [1/mV] * V) / (1 + exp(-0.1378 [1/mV] * (V + 40.14 [mV]))))
        in [1/ms]
    tau = 1 / (a + b)
        in [ms]
    inf = 1 / (1 + exp((V + 71.55 [mV]) / 7.43 [mV]))^2

#
# Background sodium current
#
# Unchanged from [2].
#
[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
#
# Rescaled, but otherwise unchanged from [2].
#
[inak]
use membrane.V, phys.FRT
use ion.K_o, sodium.Na_jn, sodium.Na_sl
INaK = INaK_jn + INaK_sl
    in [A/F]
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 = 0.99 [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]
    in [mM]

#
# Rapidly activating potassium current
#
# Unchanged from [2].
#
[ikr]
use membrane.V
IKr = gKr * sqrt(ion.K_o / 5.4 [mM]) * x * r * (V - rev.EK)
    in [A/F]
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
#
# Unchanged from [2].
#
[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 = 0.0035 [mS/uF]
    in [mS/uF]
dot(x) = (inf - x) / tau
    inf = 1 / (1 + exp(-(V + 3.8 [mV]) / 14.25 [mV]))
    tau = 990.1 [ms] / (1 + exp(-(V + 2.436 [mV]) / 14.12 [mV]))
        in [ms]

#
# Plateau potassium current
#
# No longer divided into (unused) jn and sl component, so essentially
# unchanged from [2].
#
[ikp]
use membrane.V
IKp = 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]))

#
# Transient outward potassium current (fast and slow components)
#
# Unchanged from [2].
#
[ito]
use membrane.V
Ito = Itos + Itof
    in [A/F]
Itof = gtof * xf * yf * (V - rev.EK)
    in [A/F]
Itos = gtos * xs * ys * (V - rev.EK)
    in [A/F]
gtof = if(mode.epi == 1, 0.1144 [mS/uF], 0.001404 [mS/uF])
    in [mS/uF]
gtos = if(mode.epi == 1, 0.0156 [mS/uF], 0.037596 [mS/uF])
    in [mS/uF]
x_inf = 1 / (1 + exp((V - 19 [mV]) / -13 [mV]))
y_inf = 1 / (1 + exp((V + 19.5 [mV]) / 5 [mV]))
dot(xf) = (x_inf - xf) / tau
    tau = 8.5 [ms] * exp(-((V + 45 [mV]) / 50 [mV])^2) + 0.5 [ms]
        in [ms]
dot(xs) = (x_inf - xs) / tau
    tau = 9 [ms] / (1 + exp((V + 3 [mV]) / 15 [mV])) + 0.5 [ms]
        in [ms]
dot(yf) = (y_inf - yf) / tau
    tau = 85 [ms] * exp(-(V + 40 [mV])^2 / 220 [mV^2]) + 7 [ms]
        in [ms]
dot(ys) = (y_inf - ys) / tau
    tau = 800 [ms] / (1 + exp((V + 60 [mV]) / 10 [mV])) + 30 [ms]
        in [ms]

#
# Inward recitifer potassium current
#
# Rescaled and with altered rectification.
#
[ik1]
use membrane.V, rev.EK
IK1 = gK1 * sqrt(ion.K_o / 5.4 [mM]) * inf * (V - EK)
    in [A/F]
gK1 = 0.57153 [mS/uF]
    in [mS/uF]
inf = a / (a + b)
    a = 4.0938 / (1 + exp(0.12165 [1/mV] * (V - EK - 49.9344 [mV])))
    b = (15.7197 * exp(0.06739 [1/mV] * (V - EK - 3.2571 [mV])) + exp(0.06175 [1/mV] * (V - EK - 594.31 [mV]))) / (1 + exp(-0.16285 [1/mV] * (V - EK + 14.2067 [mV])))

#
# Calcium activated chloride current
#
# Unchanged from [2].
#
[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.054813 [mS/uF]
    in [mS/uF]
KdClCa = 0.1 [mM]
    in [mM]

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

#
# L-type calcium current
#
# Updated voltage-dependent gates, and added new f2 gate.
#
[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
ICaL = ICaL_Ca + ICaL_Na + ICaL_K
    in [A/F]
ICaL_Ca = ICaL_Ca_jn + ICaL_Ca_sl
    in [A/F]
ICaL_Ca_jn = fjn_CaL * ibarca_jn * d * f * f2 * (1 - fCa_jn)
    in [A/F]
ICaL_Ca_sl = fsl_CaL * ibarca_sl * d * f * f2 * (1 - fCa_sl)
    in [A/F]
ICaL_Na = ICaL_Na_jn + ICaL_Na_sl
    in [A/F]
ICaL_Na_jn = fjn_CaL * ibarna_jn * d * f * f2 * (1 - fCa_jn)
    in [A/F]
ICaL_Na_sl = fsl_CaL * ibarna_sl * d * f * f2 * (1 - fCa_sl)
    in [A/F]
ICaL_K = ibark * d * f * f2 * (fjn_CaL * (1 - fCa_jn) + fsl_CaL * (1 - fCa_sl))
    in [A/F]
# Permeability
PCa = 1.9887e-4 [L/F/ms]
    in [L/F/ms]
PNa = 3.0375e-9 [L/F/ms]
    in [L/F/ms]
PK = 5.4675e-8 [L/F/ms]
    in [L/F/ms]
# Driving terms (gammas and 4s have all been absorbed into the P parameters)
ibarca_jn = PCa * V * F * FRT * (Ca_jn * exp(2 * V * FRT) - Ca_o) / (exp(2 * V * FRT) - 1)
    in [A/F]
ibarca_sl = PCa * V * F * FRT * (Ca_sl * exp(2 * V * FRT) - Ca_o) / (exp(2 * V * FRT) - 1)
    in [A/F]
ibarna_jn = PNa * V * F * FRT * (Na_jn * exp(V * FRT) - Na_o) / (exp(V * FRT) - 1)
    in [A/F]
ibarna_sl = PNa * V * F * FRT * (Na_sl * exp(V * FRT) - Na_o) / (exp(V * FRT) - 1)
    in [A/F]
ibark = PK * V * F * FRT * (K_i * exp(V * FRT) - K_o) / (exp(V * FRT) - 1)
    in [A/F]
# Activation
dot(d) = (inf - d) / tau
    inf = 1 / (1 + exp(-(V + 5 [mV]) / 6 [mV]))
    tau = alpha * beta + gamma
        in [ms]
    alpha = 1.4 [ms] / (1 + exp((-35 [mV] - V) / 13 [mV])) + 0.25 [ms]
        in [ms]
    beta = 1.4 / (1 + exp((V + 5 [mV]) / 5 [mV]))
    gamma = 1 [ms] / (1 + exp((50 [mV] - V) / 20 [mV]))
        in [ms]
# Voltage-dependent inactivation
dot(f) = (inf - f) / tau
    inf = 1 / (1 + exp((V + 20 [mV]) / 7 [mV]))
    tau = alpha + beta + gamma
        in [ms]
    alpha = 1102.5 [ms] * exp(-((V + 27 [mV]) / 15 [mV])^2)
        in [ms]
    beta = 200 [ms] / (1 + exp((13 [mV] - V) / 10 [mV]))
        in [ms]
    gamma = 180 [ms] / (1 + exp((V + 30 [mV]) / 10 [mV])) + 20 [ms]
        in [ms]
dot(f2) = (inf - f2) / tau
    inf = 0.67 / (1 + exp((V + 35 [mV]) / 7 [mV])) + 0.33
    tau = alpha + beta + gamma
        in [ms]
    alpha = 300 [ms] * exp(-(V + 25 [mV])^2 / 170 [mV^2])
        in [ms]
    beta = 31 [ms] / (1 + exp((25 [mV] - V) / 10 [mV]))
        in [ms]
    gamma = 16 [ms] / (1 + exp((V + 30 [mV]) / 10 [mV]))
        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
#
# Unchanged from [2].
#
[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]
INaCa_jn = cell.fjn * IbarNaCa * 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 * 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 = 4.5 [A/F]
    in [A/F]
Kdact = 0.15e-3 [mM]
    in [mM]
KmNai = 12.29 [mM]
    in [mM]
KmNao = 87.5 [mM]
    in [mM]
KmCai = 0.00359 [mM]
    in [mM]
KmCao = 1.3 [mM]
    in [mM]
ksat = 0.32
nu = 0.27

#
# Sarcolemmal calcium pump current
#
# Unchanged from [2].
#
[ipca]
IpCa = IpCa_jn + IpCa_sl
    in [A/F]
IpCa_jn = cell.fjn * Ibar_PMCA * a / (a + b)
    in [A/F]
    a = (calcium.Ca_jn / 1 [mM])^1.6
IpCa_sl = cell.fsl * Ibar_PMCA * 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]
Ibar_PMCA = 0.0673 [A/F]
    in [A/F]

#
# Background calcium current
#
# Unchanged from [2].
#
[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 = 5.513e-4 [mS/uF]
    in [mS/uF]

#
# SR calcium release
#
# Unchanged from [2].
#
[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]
    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
#
# Unchanged from [2] (except Q factor removed)
#
[serca]
use calcium.Ca_i, cabsr.Ca_sr
J_serca = 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]
    in [mM]
h = 1.787

#
# SR calcium leak
#
# Unchanged from [2].
#
[caleak]
use cabsr.Ca_sr, calcium.Ca_jn
J_SRleak = 5.348e-6 [1/ms] * (Ca_sr - Ca_jn)
    in [mM/ms]

#
# Cytosolic calcium buffering
#
# Unchanged from [2].
#
[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
J_CaB_cytosol = dot(TnCL) + dot(TnCHCa) + dot(TnCHMg) + dot(CaM) + dot(MyoCa) + dot(MyoMg) + dot(SRB)
    in [mM/ms]
# TnCl
Bmax_TnClow = 0.07 [mM]
    in [mM]
koff_tncl = 19.6e-3 [1/ms]
    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 = 0.0171 [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
#
# Unchanged from [2].
#
[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-4 [mM] * Vmyo / Vjn
    in [mM]
Bmax_SLlowsl = 37.4e-3 [mM] * Vmyo / Vsl
    in [mM]
Bmax_SLhighjn = 1.65e-4 [mM] * Vmyo / Vjn
    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
#
# Unchanged from [2].
#
[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
#
# Unchanged from [2].
#
[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 = 8.2413e-13 [L/ms]
    in [L/ms]
J_Ca_slmyo = 3.7243e-12 [L/ms]
    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]
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
#
# Unchanged from [2].
#
[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
#
# Different Js compared to [2]
#
[sodium]
use cell.C, phys.F
use cell.Vjn, cell.Vsl, cell.Vmyo
J_Na_jnsl = 1.8313e-14 [m^3/s]
    in [m^3/s]
J_Na_slmyo = 1.6386e-12 [m^3/s]
    in [m^3/s]
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]

#
# Potassium concentration
#
# Unchanged from [2].
#
[potassium]
I_K_tot = ito.Ito + ikr.IKr + iks.IKs + ik1.IK1 - 2 * inak.INaK + ical.ICaL_K + ikp.IKp
    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()