stewart-2009.mmt

[[model]]
name: stewart-2009
version: 20240904
mmt_authors: Michael Clerx
display_name: Stewart et al., 2009
desc: """
    The 2009 model of the human Purinkje action potential by Stewart et
    al. [1], based on the ventricular model by Ten Tusscher et al. [2].

    This implementation was created from the CellML code made available with
    the original publication [2]. Several "cosmetic" changes were made without
    affecting the equations, which were tested against the original CellML and
    which matched to within rounding error precision.
    - The state variables were ordered as in the publication
    - The components were ordered as in the publication, and were annotated
      with links to the relevant document section.
    In addition, incorrect units in F, R, Cm, gCaL, Vc, Vss, and Vsr were
    corrected.

    The stimulus was set to 0.5ms and approximately twice the threshold for
    (immediate) depolarisation.

    [1] Stewart, P., Aslanidi, O. V., Noble, D., Noble, P. J., Boyett, M. R., &
        Zhang, H. (2009). Mathematical model of the electrical action potential
        of Purkinje fibre cells. Philosophical Transactions of the Royal
        Society, 367(1896), 2225-2255.
        https://doi.org/10.1098/rsta.2008.0283

    [2] Ten Tusscher, K. H. W. J., & Panfilov, A. V. (2006). Alternans and
        spiral breakup in a human ventricular tissue model. American Journal of
        Physiology. Heart and Circulatory Physiology, 291(3), H1088-H1100.
        https://doi.org/10.1152/ajpheart.00109.2006

    [3] https://models.physiomeproject.org/exposure/38cf8387b0707f0ef6947f009710aeb5
        Retrieved on 2021-10-27

    An extract from the CellML file's meta data follows below:
    ---------------------------------------------------------------------------
    Penny Noble
    DPAG, University of Oxford

    Model Status
    This CellML ersion of the model has been checked in COR and OpenCell. The
    units are consistent and the model runs to recreate the published results.

    Model Structure
    In the paper described here, Philip Stewart and colleagues present a
    details of their newly developed model for the human Purkinje cell
    including validation against experimental data. Ionic mechanisms underlying
    the heterogeneity between the Purkinje fibre and ventricular action
    potentials in humans and other species were analysed. The newly developed
    Purkinje fibre cell model adds a new member to the family of human cardiac
    cell models developed previously for the sino-atrial node, atrial and
    ventricular cells, which can be incorporated into an anatomical model of
    human heart with details of its electrophysiological heterogeneity and
    anatomical complexity.
    """
# Initial values
membrane.V    = -7.12864384994752527e+01
calcium.Ca_i  =  1.02363913704157998e-04
calcium.Ca_SR =  3.14149868138687083e+00
calcium.Ca_ss =  3.81250245527617196e-04
sodium.Na_i   =  8.80505373054986329e+00
potassium.K_i =  1.36773426842998674e+02
ina.m         =  2.83293473203283658e-02
ina.h         =  2.51110676365191021e-01
ina.j         =  2.65680422809985550e-01
ikr.Xr1       =  1.34628009070351228e-02
ikr.Xr2       =  3.32697379577155727e-01
iks.Xs        =  8.40127007165385360e-03
ito.r         =  8.84290083740064937e-04
ito.s         =  9.68749754333663282e-01
ical.d        =  2.16290842830038621e-04
ical.f1       =  9.68969771502977917e-01
ical.f2       =  9.96328488888527430e-01
ical.fCass    =  9.99935288766206076e-01
if.y          =  3.84052487230444398e-02
irel.R_prime  =  9.84086493577585375e-01

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

#
# Membrane potential
#
[membrane]
use stimulus.i_stim
dot(V) = -(i_ion + i_stim + i_diff)
    in [mV]
    label membrane_potential
i_ion = sodium.INa_tot + potassium.IK_tot + calcium.ICa_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]
amplitude = -10 [A/F] * 2
    in [A/F]

#
# (a) Inward rectifier current, IK1
# Appendix page 2242
#
[ik1]
use membrane.V
inf = 1 / (1 + exp(0.1 [1/mV] * (V + 75.44 [mV])))
g_K1 = 0.065 [mS/uF]
    in [mS/uF]
i_K1 = g_K1 * inf * (V - 8 [mV] - nernst.E_K)
    in [A/F]

#
# (b) Transient outward current, ITo
# Appendix page 2242
#
[ito]
use membrane.V
dot(r) = (r_inf - r) / tau_r
    r_inf = 1 / (1 + exp((20 [mV] - V) / 13 [mV]))
    tau_r = 10.45 [ms] * exp(-(V + 40 [mV])^2 / 1800 [mV^2]) + 7.3 [ms]
        in [ms]
dot(s) = (s_inf - s) / tau_s
    s_inf = 1 / (1 + exp((V + 27 [mV]) / 13 [mV]))
    tau_s = 85 [ms] * exp(-(V + 25 [mV])^2 / 320 [mV^2]) + 5 [ms] / (1 + exp((V - 40 [mV]) / 5 [mV])) + 42 [ms]
        in [ms]
g_to = 0.08184 [mS/uF]
    in [mS/uF]
i_to = g_to * r * s * (V - nernst.E_K)
    in [A/F]

#
# (c) Sustained current, Isus
# Appendix page 2243
#
[isus]
use membrane.V
a = 1 / (1 + exp((5 [mV] - V) / 17 [mV]))
g_sus = 0.0227 [mS/uF]
    in [mS/uF]
i_sus = g_sus * a * (V - nernst.E_K)
    in [A/F]

#
# (d) Hyperpolarization-activated current, If
# Appendix page 2243
# Described in detail on page 2229
#
[if]
use membrane.V
dot(y) = (y_inf - y) / tau_y
    alpha_y = 1 [1/ms] * exp(-2.9 - 0.04 [1/mV] * V)
        in [1/ms]
    beta_y = 1 [1/ms] * exp(3.6 + 0.11 [1/mV] * V)
        in [1/ms]
    tau_y = 4000 / (alpha_y + beta_y)
        in [ms]
    y_inf = 1 / (1 + exp((V + 80.6 [mV]) / 6.8 [mV]))
g_f_K = 0.0234346 [mS/uF]
    in [mS/uF]
g_f_Na = 0.0145654 [mS/uF]
    in [mS/uF]
i_f_K = y * g_f_K * (V - nernst.E_K)
    in [A/F]
i_f_Na = y * g_f_Na * (V - nernst.E_Na)
    in [A/F]
i_f = i_f_Na + i_f_K
    in [A/F]

#
# (e) Fast sodium current, INa
# Appendix page 2243
#
[ina]
use membrane.V
dot(m) = (inf - m) / tau
    alpha = 1 / (1 + exp((-60 [mV] - V) / 5 [mV]))
    beta = 0.1 / (1 + exp((V + 35 [mV]) / 5 [mV])) + 0.1 / (1 + exp((V - 50 [mV]) / 200 [mV]))
    tau = 1 [ms] * alpha * beta
        in [ms]
    inf = 1 / (1 + exp((-56.86 [mV] - V) / 9.03 [mV]))^2
dot(h) = (inf - h) / tau
    alpha = if(V < -40 [mV], 0.057 [1/ms] * exp(-(V + 80 [mV]) / 6.8 [mV]), 0 [1/ms])
        in [1/ms]
    beta = if(V < -40 [mV],
              2.7 [1/ms] * exp(0.079 [1/mV] * V) + 310000 [1/ms] * exp(0.3485 [1/mV] * V),
              0.77 [1/ms] / (0.13 * (1 + exp((V + 10.66 [mV]) / -11.1 [mV]))))
        in [1/ms]
    tau = 1 / (alpha + beta)
        in [ms]
    inf = 1 / (1 + exp((V + 71.55 [mV]) / 7.43 [mV]))^2
dot(j) = (inf - j) / tau
    alpha = if(V < -40 [mV],
               (-25428 [1/ms] * exp(0.2444 [1/mV] * V) - 6.948e-6 [1/ms] * exp(-0.04391 [1/mV] * V)) * (V + 37.78 [mV]) / 1 [mV] / (1 + exp(0.311 [1/mV] * (V + 79.23 [mV]))),
               0 [1/ms])
        in [1/ms]
    beta = if(V < -40 [mV],
              0.02424 [1/ms] * exp(-0.01052 [1/mV] * V) / (1 + exp(-0.1378 [1/mV] * (V + 40.14 [mV]))),
              0.6 [1/ms] * exp(0.057 [1/mV] * V) / (1 + exp(-0.1 [1/mV] * (V + 32 [mV]))))
        in [1/ms]
    tau = 1 / (alpha + beta)
        in [ms]
    inf = 1 / (1 + exp((V + 71.55 [mV]) / 7.43 [mV]))^2
g_Na = 130.5744 [mS/uF]
    in [mS/uF]
i_Na = g_Na * m^3 * h * j * (V - nernst.E_Na)
    in [A/F]

#
# (f) L-type calcium current, ICaL
# Appendix page 2244
#
[ical]
use membrane.V
use phys.FRT, phys.FFRT
use calcium.Ca_ss, extra.Ca_o
dot(d) = (d_inf - d) / tau_d
    alpha_d = 1.4 / (1 + exp((-35 [mV] - V) / 13 [mV])) + 0.25
    beta_d = 1.4 / (1 + exp((V + 5 [mV]) / 5 [mV]))
    gamma_d = 1 [ms] / (1 + exp((50 [mV] - V) / 20 [mV]))
        in [ms]
    tau_d = 1 [ms] * alpha_d * beta_d + gamma_d
        in [ms]
    d_inf = 1 / (1 + exp((-8 [mV] - V) / 7.5 [mV]))
dot(f1) = (f_inf - f1) / tau_f
    f_inf = 1 / (1 + exp((V + 20 [mV]) / 7 [mV]))
    tau_f = 1102.5 [ms] * exp(-(V + 27 [mV])^2 / 225 [mV^2]) + 200 [ms] / (1 + exp((13 [mV] - V) / 10 [mV])) + 180 [ms] / (1 + exp((V + 30 [mV]) / 10 [mV])) + 20 [ms]
        in [ms]
dot(f2) = (f2_inf - f2) / tau_f2
    f2_inf = 0.67 / (1 + exp((V + 35 [mV]) / 7 [mV])) + 0.33
    tau_f2 = 562 [ms] * exp(-(V + 27 [mV])^2 / 240 [mV^2]) + 31 [ms] / (1 + exp((25 [mV] - V) / 10 [mV])) + 80 [ms] / (1 + exp((V + 30 [mV]) / 10 [mV]))
        in [ms]
dot(fCass) = (fCass_inf - fCass) / tau_fCass
    fCass_inf = 0.6 / (1 + (calcium.Ca_ss / 0.05 [mM])^2) + 0.4
    tau_fCass = 80 [ms] / (1 + (calcium.Ca_ss / 0.05 [mM])^2) + 2 [ms]
        in [ms]
g_CaL = 0.0398 [L/F/s]
    in [L/F/s]
i_CaL = if(abs(V - 15 [mV]) < 1e-6 [mV], a * (b - Ca_o) / (2 * FRT), numer / denom)
    in [A/F]
    a = g_CaL * d * f1 * f2 * fCass * 4 * FFRT
        in [1/mM/ms]
    b = 0.25 * Ca_ss
        in [mM]
    numer = a * (V - 15 [mV]) * (b * exp((V - 15 [mV]) * 2 * FRT) - Ca_o)
        in [A/F]
    denom = (exp(2 * (V - 15 [mV]) * FRT) - 1)

#
# (g) Slow delayed rectifier current, IKs
# Appendix page 2245
# Also called the "slow time-dependent potassium current"
#
[iks]
use membrane.V
dot(Xs) = (inf - Xs) / tau
    alpha = 1400 / sqrt(1 + exp((5 [mV] - V) / 6 [mV]))
    beta = 1 / (1 + exp((V - 35 [mV]) / 15 [mV]))
    tau = 1 [ms] * alpha * beta + 80 [ms]
        in [ms]
    inf = 1 / (1 + exp((-5 [mV] - V) / 14 [mV]))
g_Ks = 0.2352 [mS/uF]
    in [mS/uF]
i_Ks = g_Ks * Xs^2 * (V - nernst.E_Ks)
    in [A/F]

#
# (h) Rapid delayed rectifier current, IKr
# Appendix page 2246
# Also called the "rapid time-dependent potassium current"
#
[ikr]
use membrane.V
dot(Xr1) = (xr1_inf - Xr1) / tau_xr1
    alpha_xr1 = 450 / (1 + exp((-45 [mV] - V) / 10 [mV]))
    beta_xr1 = 6 / (1 + exp((V + 30 [mV]) / 11.5 [mV]))
    tau_xr1 = 1 [ms] * alpha_xr1 * beta_xr1
        in [ms]
    xr1_inf = 1 / (1 + exp((-26 [mV] - V) / 7 [mV]))
dot(Xr2) = (xr2_inf - Xr2) / tau_xr2
    alpha_xr2 = 3 / (1 + exp((-60 [mV] - V) / 20 [mV]))
    beta_xr2 = 1.12 / (1 + exp((V - 60 [mV]) / 20 [mV]))
    tau_xr2 = 1 [ms] * alpha_xr2 * beta_xr2
        in [ms]
    xr2_inf = 1 / (1 + exp((V + 88 [mV]) / 24 [mV]))
g_Kr = 0.0918 [mS/uF]
    in [mS/uF]
    label g_Kr
i_Kr = g_Kr * sqrt(extra.K_o / 5.4 [mM]) * Xr1 * Xr2 * (V - nernst.E_K)
    in [A/F]

#
# (i) Na/Ca exchange current, INaCa
# Appendix page 2246
#
[inaca]
use membrane.V, phys.FRT
use extra.Na_o, extra.Ca_o
use sodium.Na_i, calcium.Ca_i
alpha = 2.5
gamma = 0.35
K_sat = 0.1
K_NaCa = 1000 [A/F]
    in [A/F]
Km_Nai = 87.5 [mM]
    in [mM]
Km_Ca = 1.38 [mM]
    in [mM]
i_NaCa = K_NaCa * (exp(gamma * V * FRT) * Na_i^3 * Ca_o - exp((gamma - 1) * V * FRT) * Na_o^3 * Ca_i * alpha) / ((Na_o^3 + Km_Nai^3) * (Km_Ca + Ca_o) * (1 + K_sat * exp((gamma - 1) * V * FRT)))
    in [A/F]

#
# (j) Na/K pump current, INaK
# Appendix page 2247
#
[inak]
use membrane.V, phys.FRT
use extra.K_o, sodium.Na_i
K_mNa = 40 [mM]
    in [mM]
K_mk = 1 [mM]
    in [mM]
P_NaK = 2.724 [A/F]
    in [A/F]
i_NaK = P_NaK * K_o / (K_o + K_mk) * Na_i / (Na_i + K_mNa) / (1 + 0.1245 * exp(-0.1 * V * FRT) + 0.0353 * exp(-V * FRT))
    in [A/F]

#
# (j) Calcium pump current, IpCa
# Appendix page 2247
#
[ipca]
K_pCa = 0.0005 [mM]
    in [mM]
g_pCa = 0.1238 [A/F]
    in [A/F]
i_p_Ca = g_pCa * calcium.Ca_i / (calcium.Ca_i + K_pCa)
    in [A/F]

#
# (j) Potassium pump current, IpK
# Appendix page 2247
#
[ipk]
use membrane.V
g_pK = 0.0146 [mS/uF]
    in [mS/uF]
i_p_K = g_pK * (V - nernst.E_K) / (1 + exp((25 [mV] - V) / 5.98 [mV]))
    in [A/F]

#
# (k) Sodium background current, IbNa
# Appendix page 2247
#
[ibna]
g_bna = 0.00029 [mS/uF]
    in [mS/uF]
i_b_Na = g_bna * (membrane.V - nernst.E_Na)
    in [A/F]

#
# (k) Calcium background current, IbCa
# Appendix page 2247
#
[ibca]
g_bca = 0.000592 [mS/uF]
    in [mS/uF]
i_b_Ca = g_bca * (membrane.V - nernst.E_Ca)
    in [A/F]

#
# (l) Calcium uptake into the SR (SERCA)
#
[iup]
use calcium.Ca_i
Vmax_up = 0.006375 [mM/ms]
    in [mM/ms]
K_up = 0.00025 [mM]
    in [mM]
i_up = Vmax_up / (1 + K_up^2 / Ca_i^2)
    in [mM/ms]

#
# (l) Calcium release from the SR (RyR)
#
[irel]
use calcium.Ca_ss, calcium.Ca_SR
k1_prime = 0.15 [1/mM^2/ms]
    in [1/mM^2/ms]
k2_prime = 0.045 [1/mM/ms]
    in [1/mM/ms]
kcasr = max_sr - (max_sr - min_sr) / (1 + (EC / Ca_SR)^2)
    max_sr = 2.5
    min_sr = 1
    EC = 1.5 [mM]
        in [mM]
k1 = k1_prime / kcasr
    in [1/mM^2/ms]
k2 = k2_prime * kcasr
    in [1/mM/ms]
k3 = 0.06 [1/ms]
    in [1/ms]
k4 = 0.005 [1/ms]
    in [1/ms]
dot(R_prime) = -k2 * Ca_ss * R_prime + k4 * (1 - R_prime)
O = k1 * Ca_ss^2 * R_prime / (k3 + k1 * Ca_ss^2)
V_rel = 0.102 [mS/uF]
    in [mS/uF]
i_rel = V_rel * O * (Ca_SR - Ca_ss)
    in [mM/ms]

#
# (l) Calcium leak from SR
#
[ileak]
use calcium.Ca_SR, calcium.Ca_i
V_leak = 0.00036 [mS/uF]
    in [mS/uF]
i_leak = V_leak * (Ca_SR - Ca_i)
    in [mM/ms]

#
# (l) Diffusion from SS to Bulk
#
[ixfer]
use calcium.Ca_i, calcium.Ca_ss
V_xfer = 0.0038 [mS/uF]
    in [mS/uF]
i_xfer = V_xfer * (Ca_ss - Ca_i)
    in [mM/ms]

#
# (l) Calcium dynamics
# Appendix page 2247
#
[calcium]
use membrane.V
use cell.V_c, cell.V_sr, cell.V_ss, cell.Cm
use ileak.i_leak, iup.i_up, irel.i_rel, ixfer.i_xfer
# Buffering constants
Buf_c = 0.2 [mM]
    in [mM]
Buf_ss = 0.4 [mM]
    in [mM]
Buf_sr = 10 [mM]
    in [mM]
K_buf_c = 0.001 [mM]
    in [mM]
K_buf_ss = 0.00025 [mM]
    in [mM]
K_buf_sr = 0.3 [mM]
    in [mM]
# Buffered concentrations
Ca_i_bufc = 1 / (1 + Buf_c * K_buf_c / (Ca_i + K_buf_c)^2)
Ca_ss_bufss = 1 / (1 + Buf_ss * K_buf_ss / (Ca_ss + K_buf_ss)^2)
Ca_sr_bufsr = 1 / (1 + Buf_sr * K_buf_sr / (Ca_SR + K_buf_sr)^2)
# Free calcium in cytosol, SS, and SR
dot(Ca_i) = Ca_i_bufc * ((i_leak - i_up) * V_sr / V_c + i_xfer - ICa_cyt * Cm / (2 * V_c * phys.F))
    in [mM]
dot(Ca_ss) = Ca_ss_bufss * (-ical.i_CaL * Cm / (2 * V_ss * phys.F) + i_rel * V_sr / V_ss - i_xfer * V_c / V_ss)
    in [mM]
dot(Ca_SR) = Ca_sr_bufsr * (i_up - (i_rel + i_leak))
    in [mM]
# Total calcium currents
ICa_cyt = ibca.i_b_Ca + ipca.i_p_Ca - 2 * inaca.i_NaCa
    in [A/F]
ICa_tot = ICa_cyt + ical.i_CaL
    in [A/F]

#
# (m) Sodium dynamics
# Appendix page 2248
#
[sodium]
INa_tot = ina.i_Na + ibna.i_b_Na + if.i_f_Na + 3 * inak.i_NaK + 3 * inaca.i_NaCa
    in [A/F]
dot(Na_i) = -INa_tot / (cell.V_c * phys.F) * cell.Cm
    in [mM]

#
# (m) Potassium dynamics
# Appendix page 2248
#
[potassium]
IK_tot = ik1.i_K1 + ito.i_to + if.i_f_K + isus.i_sus + ikr.i_Kr + iks.i_Ks + ipk.i_p_K - 2 * inak.i_NaK
    in [A/F]
dot(K_i) = -IK_tot / (cell.V_c * phys.F) * cell.Cm  # cell.Cm is not in the appendix!
    in [mM]

#
# Physical constants
#
[phys]
R = 8.314472 [J/mol/K]
    in [J/mol/K]
T = 310 [K]
    in [K]
F = 9.64853414999999950e+01 [C/mmol]
    in [C/mmol]
RTF = R * T / F
    in [mV]
FRT = F / R / T
    in [1/mV]
FFRT = F * FRT
    in [C/mmol/mV]

#
# Cell geometry
#
[cell]
Cm = 185 [pF]
    in [pF]
    desc: Cell capacitance
V_c = 16404 [um^3]
    in [um^3]
    desc: Bulk cytoplasm volume
V_ss = 54.68 [um^3]
    in [um^3]
    desc: Dyadic (junctional) subspace volume
V_sr = 1094 [um^3]
    in [um^3]
    desc: Sarcoplasmic reticulum volume

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

#
# Reversal potentials
#
[nernst]
use extra.K_o, extra.Na_o, extra.Ca_o
use potassium.K_i, sodium.Na_i, calcium.Ca_i
E_Ca = 0.5 * phys.RTF * log(Ca_o / Ca_i)
    in [mV]
E_K = phys.RTF * log(K_o / K_i)
    in [mV]
E_Ks = phys.RTF * log((K_o + P_kna * Na_o) / (K_i + P_kna * Na_i))
    in [mV]
E_Na = phys.RTF * log(Na_o / Na_i)
    in [mV]
P_kna = 0.03

[[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)

# Run simulation
d = s.run(3000)

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