bai-2018.mmt

[[model]]
name: bai-2018
version: 20240904
mmt_authors: Michael Clerx
display_name: Bai et al., 2018
desc: """
    The human left-atrial AP model by Bai et al. [1], based on the
    ventricular model by Ten Tusscher and Panfilov [2].

    This implementation is based on the author-provided CellML [3], but
    includes the unit fixes described in [4]. The implementation was verified
    numerically by comparing the calculates derivatives to the values from the
    author-provided CellML. Differences from [2] are indicated throughout.

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

    References:

    [1] Bai, J. Gladding, P. A., Stiles, M. K., Fedorov, V. V. & Zhao, J.
        (2018) Ionic and cellular mechanisms underlying TBX5/PITX2
        insufficiency-induced atrial fibrillation: Insights from mathematical
        models of human atrial cells. Scientific Reports 8:15642
        https://doi.org/10.1038/s41598-018-33958-y

    [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.cellml.org/workspace/520

    [4] Clerx, M., Barral, Y., Agrawal, A. & Mirams, G. R. (2020)
        TNNP unit fixes
        https://models.cellml.org/workspace/604/file/de5a4e600b57c8b9b0bd477695648c82351bc6dd/TNNP_unit_fixes.pdf

    """
# Initial values
membrane.V   = -77.6
calcium.Cai  = 7e-5
calcium.CaSR = 1.3
calcium.CaSS = 7e-5
sodium.Nai   = 7.67
potassium.Ki = 138.3
ina.m        = 0
ina.h        = 0.75
ina.j        = 0.75
ikr.xr1      = 0
ikr.xr2      = 1
iks.xs       = 0
ito.r        = 0
ito.s        = 1
ikur.x       = 3.824458e-4
ikur.y       = 0.9597222
ical.d       = 0
ical.f       = 1
ical.f2      = 1
ical.fCaSS   = 1
jrel.R       = 1

#
# Simulator 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)
    in [mV]
    label membrane_potential
i_ion = (+ ina.INa
         + ik1.IK1
         + ikr.IKr
         + iks.IKs
         + ito.Ito
         + ikur.IKur
         + ical.ICaL
         + inak.INaK
         + inaca.INaCa
         + ipca.IpCa
         + ipk.IpK
         + icab.ICab
         + inab.INab
        )
    in [A/F]
    label cellular_current

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

#
# Cell properties
#
[cell]
Vc = 16404 [um^3]
    in [um^3]
    desc: Bulk cytoplasm volume
Vss = 54.68 [um^3]
    in [um^3]
    desc: Dyadic (junctional) subspace volume
Vsr = 1094 [um^3]
    in [um^3]
    desc: Sarcoplasmic reticulum volume
Cm = 185 [pF]
    in [pF]
    desc: Cell capacitance
    label membrane_capacitance

#
# Physical constants
#
[phys]
F = 96.485 [C/mmol]
    in [C/mmol]
R = 8.314 [J/mol/K]
    in [J/mol/K]
T = 310 [K]
    in [K]
RTF = R * T / F
    in [mV]
FRT = F / (R * T)
    in [1/mV]
FFRT = F * FRT
    in [C/mmol/mV]

#
# Fast sodium current
#
[ina]
use membrane.V, nernst.ENa
INa = gNa * m^3 * h * j * (V - ENa)
    in [A/F]
gNa = 14.838 [mS/uF]
    in [mS/uF]
dot(m) = (inf - m) / tau
    inf = 1 / (1 + exp((-56.86 [mV] - V) / 9.03 [mV]))^2
    tau = 1 [ms] * alpha * beta
        in [ms]
    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]))
dot(h) = (inf - h) / tau
    inf = 1 / (1 + exp((V + 71.55 [mV]) / 7.43 [mV]))^2
    tau = 1 / (alpha + beta)
        in [ms]
    alpha = if(V < -40 [mV],
                 0.057 [mS/uF] * exp(-(V + 80 [mV]) / 6.8 [mV]),
                 0 [mS/uF])
        in [1/ms]
    beta = if(V < -40 [mV],
                2.7 [mS/uF] * exp(0.079 [1/mV] * V) + 310000 [mS/uF] * exp(0.3485 [1/mV] * V),
                0.77 [mS/uF] / (0.13 * (1 + exp((V + 10.66 [mV]) / -11.1 [mV]))))
        in [1/ms]
dot(j) = (inf - j) / tau
    inf = 1 / (1 + exp((V + 71.55 [mV]) / 7.43 [mV]))^2
    tau = 1 / (alpha + beta)
        in [ms]
    alpha = if(V < -40 [mV],
               (-25428 [mS/uF] * exp(0.2444 [1/mV] * V) - 6.948e-6 [mS/uF] * exp(-0.04391 [1/mV] * V)) * (V + 37.78 [mV]) / 1 [mV] / (1 + exp(0.311 [1/mV] * (V + 79.23 [mV]))),
               0 [mS/uF])
        in [1/ms]
    beta = if(V < -40 [mV],
              0.02424 [mS/uF] * exp(-0.01052 [1/mV] * V) / (1 + exp(-0.1378 [1/mV] * (V + 40.14 [mV]))),
              0.6 [mS/uF] * exp(0.057 [1/mV] * V) / (1 + exp(-0.1 [1/mV] * (V + 32 [mV]))))
        in [1/ms]

#
# Inward rectifier potassium current
#
[ik1]
use membrane.V, nernst.EK
IK1 = gK1 * inf * (V - EK)
    in [A/F]
gK1 = 1.081 [mS/uF]  # Rescaled from original, and dropped EK dependence
    in [mS/uF]
inf = alpha / (alpha + beta)
    alpha = 0.1 [1/ms] / (1 + exp(0.06 [1/mV] * (V - EK - 200 [mV])))
        in [1/ms]
    beta = (3 [1/ms] * exp(0.0002 [1/mV] * (V - EK + 100 [mV])) + 1 [1/ms] * exp(0.1 [1/mV] * (V - EK - 10 [mV]))) / (1 + exp(-0.5 [1/mV] * (V - EK)))
        in [1/ms]

#
# Rapid time-dependent potassium current
#
[ikr]
use membrane.V, nernst.EK
IKr = gKr * sqrt(extra.Ko / 5.4 [mM]) * xr1 * xr2 * (V - EK)
    in [A/F]
gKr = 0.1989 [mS/uF]    # Rescaled
    in [mS/uF]
    label g_Kr
dot(xr1) = (inf - xr1) / tau
    inf = 1 / (1 + exp((-26 [mV] - V) / 7 [mV]))
    tau = 1 [ms] * alpha * beta
        in [ms]
    alpha = 450 / (1 + exp((-45 [mV] - V) / 10 [mV]))
    beta = 6 / (1 + exp((V + 30 [mV]) / 11.5 [mV]))
dot(xr2) = (inf - xr2) / tau
    inf = 1 / (1 + exp((V + 88 [mV]) / 24 [mV]))
    tau = 1 [ms] * alpha * beta
        in [ms]
    alpha = 3 / (1 + exp((-60 [mV] - V) / 20 [mV]))
    beta = 1.12 / (1 + exp((V - 60 [mV]) / 20 [mV]))

#
# Slow time-dependent potassium current
#
[iks]
use membrane.V, nernst.EKs
IKs = gKs * xs^2 * (V - EKs)
    in [A/F]
gKs = 0.784 [mS/uF]     # Rescaled
    in [mS/uF]
dot(xs) = (inf - xs) / tau
    inf = 1 / (1 + exp((-5 [mV] - V) / 14 [mV]))
    tau = 1 [ms] * alpha * beta + 80 [ms]
        in [ms]
    alpha = 1400 / sqrt(1 + exp((5 [mV] - V) / 6 [mV]))
    beta = 1 / (1 + exp((V - 35 [mV]) / 15 [mV]))

#
# Transient outward current
#
[ito]
use membrane.V, nernst.EK
Ito = gto * r * s * (V - EK)
    in [A/F]
gto = 0.588 [mS/uF]     # Rescaled
    in [mS/uF]
dot(r) = (inf - r) / tau
    inf = 1 / (1 + exp((20 [mV] - V) / 6 [mV]))
    tau = 9.5 [ms] * exp(-(V + 40 [mV])^2 / 1800 [mV^2]) + 0.8 [ms]
        in [ms]
dot(s) = (inf - s) / tau
    # The epicardial equations from [2] are used
    inf = 1 / (1 + exp((V + 20 [mV]) / 5 [mV]))
    tau = 85 [ms] * exp(-(V + 45 [mV])^2 / 320 [mV^2]) + 5 [ms] / (1 + exp((V - 20 [mV]) / 5 [mV])) + 3 [ms]
        in [ms]

#
# Ultra-rapid potassium current
#
# Equations 4 to 10 in [1].
#
[ikur]
use membrane.V, nernst.EK
IKur = gKur * x * y * (V - EK)
    in [A/F]
    label I_Kur
gKur = 0.045 [mS/uF]
    in [mS/uF]
    label g_Kur
dot(x) = (inf - x) / tau
    inf = 1 / (1 + exp((V + 6 [mV]) / -8.6 [mV]))
    tau = 9 [ms] / (1 + exp((V + 6 [mV]) / 12 [mV])) + 0.5 [ms]
        in [ms]
    # Note: the paper says exp((V + 5 [mV]
dot(y) = (inf - y) / tau
    inf = 1 / (1 + exp((V + 7.5 [mV]) / 10 [mV]))
    tau = 590 [ms] / (1 + exp((V + 60 [mV]) / 10 [mV])) + 3050 [ms]
        in [ms]

#
# L-type calcium current
#
[ical]
use membrane.V
use phys.FRT, phys.FFRT
use extra.Cao, calcium.Cai, calcium.CaSS
ICaL = gCaL * d * f * f2 * fCaSS * 4 * (V - 15 [mV]) * FFRT * (0.25 * CaSS * exp(2 * (V - 15 [mV]) * FRT) - Cao) / (exp(2 * (V - 15 [mV]) * FRT) - 1)
    in [A/F]
gCaL = 0.0597 [L/F/s]  # Rescaled
    in [L/F/s]
dot(d) = (inf - d) / tau
    inf = 1 / (1 + exp((-8 [mV] - V) / 7.5 [mV]))
    tau = 1 [ms] * (alpha * beta + gamma)
        in [ms]
    alpha = 1.4 / (1 + exp((-35 [mV] - V) / 13 [mV])) + 0.25
    beta = 1.4 / (1 + exp((V + 5 [mV]) / 5 [mV]))
    gamma = 1 / (1 + exp((50 [mV] - V) / 20 [mV]))
dot(f) = (inf - f) / tau
    inf = 1 / (1 + exp((V + 20 [mV]) / 7 [mV]))
    tau = 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) = (inf - f2) / tau
    inf = 0.67 / (1 + exp((V + 35 [mV]) / 7 [mV])) + 0.33
    # The tau equation has been reparameterised
    tau = 600 [ms] * exp(-(V + 25 [mV])^2 / 170 [mV^2]) + 31 [ms] / (1 + exp((25 [mV] - V) / 10 [mV])) + 16 [ms] / (1 + exp((V + 30 [mV]) / 10 [mV]))
        in [ms]
dot(fCaSS) = (inf - fCaSS) / tau
    inf = 0.6 / (1 + (CaSS / 0.05 [mM])^2) + 0.4
    tau = 80 [ms] / (1 + (CaSS / 0.05 [mM])^2) + 2 [ms]
        in [ms]

#
# Sodium-potassium pump
#
[inak]
use membrane.V, phys.FRT
use extra.Ko, sodium.Nai
INaK = PNaK * Ko / (Ko + K_mk) * Nai / (Nai + K_mNa) / (1 + 0.1245 * exp(-0.1 * V * FRT) + 0.0353 * exp(-V * FRT))
    in [A/F]
PNaK = 1.9068 [A/F]  # Rescaled
    in [A/F]
K_mNa = 40 [mM]
    in [mM]
K_mk = 1 [mM]
    in [mM]

#
# Sodium-calcium exchanger
#
[inaca]
use membrane.V, phys.FRT
use extra.Nao, extra.Cao
use sodium.Nai, calcium.Cai
INaCa = (K_NaCa
         * (exp(gamma * V * FRT) * Nai^3 * Cao - exp((gamma - 1) * V * FRT) * Nao^3 * Cai * alpha)
         / ((Km_Nai^3 + Nao^3) * (Km_Ca + Cao) * (1 + K_sat * exp((gamma - 1) * V * FRT))))
    in [A/F]
K_NaCa = 700 [A/F]  # Rescaled
    in [A/F]
Km_Ca = 1.38 [mM]
    in [mM]
Km_Nai = 87.5 [mM]
    in [mM]
K_sat = 0.08    # Rescaled
alpha = 2.5
gamma = 0.35

#
# Calcium pump
#
[ipca]
use calcium.Cai
IpCa = gpCa * Cai / (Cai + KpCa)
    in [A/F]
gpCa = 0.1238 [A/F]
    in [A/F]
KpCa = 0.0005 [mM]
    in [mM]

#
# Potassium pump
#
[ipk]
use membrane.V, nernst.EK
IpK = gpK * (V - EK) / (1 + exp((25 [mV] - V) / 5.98 [mV]))
    in [A/F]
gpK = 0.0146 [mS/uF]
    in [mS/uF]

#
# Background calcium current
#
[icab]
use membrane.V, nernst.ECa
ICab = gCab * (V - ECa)
    in [A/F]
gCab = 0.000592 [mS/uF]
    in [mS/uF]

#
# Background sodium current
#
[inab]
use membrane.V, nernst.ENa
INab = gNab * (V - ENa)
    in [A/F]
gNab = 0.00029 [mS/uF]
    in [mS/uF]

#
# External concentrations
#
[extra]
Cao = 2 [mM]
    in [mM]
Nao = 140 [mM]
    in [mM]
Ko = 5.4 [mM]
    in [mM]

#
# Nernst/Reversal potentials
#
[nernst]
use extra.Cao, calcium.Cai
use extra.Nao, sodium.Nai
use extra.Ko, potassium.Ki
use phys.RTF
ECa = RTF * log(Cao / Cai) * 0.5
    in [mV]
ENa = RTF * log(Nao / Nai)
    in [mV]
EK = RTF * log(Ko / Ki)
    in [mV]
EKs = RTF * log((Ko + P_kna * Nao) / (Ki + P_kna * Nai))
    in [mV]
P_kna = 0.03

#
# Calcium release from the SR (RyR and leak)
# Updated by equations 2 and 3 in [1]
#
[jrel]
use membrane.V
use calcium.CaSS, calcium.CaSR
Jrel = Jcicr + Jleak
    in [mM/ms]
Jcicr = Vrel * O * (CaSR - CaSS)
    in [mM/ms]
Jleak = Vleak * R * (CaSR - CaSS)
    in [mM/ms]
Vrel = 0.102 [1/ms]
    in [1/ms]
Vleak = 0.00036 [1/ms]
    in [1/ms]
kcasr = max_sr - (max_sr - min_sr) / (1 + (EC / CaSR)^2)
    max_sr = 2.5
    min_sr = 1
    EC = 1.5 [mM]
        in [mM]
k1 = 0.15 [1/mM^2/ms] / kcasr
    in [1/mM^2/ms]
k2 = 0.045 [1/mM/ms] * kcasr
    in [1/mM/ms]
k3 = 0.06 [mS/uF]
    in [mS/uF]
k4 = 0.005 [mS/uF]
    in [mS/uF]
dot(R) = -k2 * CaSS * R + k4 * (1 - R)
O = k1 * CaSS^2 * R / (k3 + k1 * CaSS^2)

#
# Calcium uptake into the SR (SERCA)
#
[jup]
use calcium.Cai
Jup = Vmax_up / (1 + K_up^2 / Cai^2)
    in [mM/ms]
Vmax_up = 0.006375 [mM/ms]
    in [mM/ms]
K_up = 0.0005 [mM]      # Doubled
    in [mM]

#
# Diffusion from SS to Bulk
#
[jxfer]
use calcium.Cai, calcium.CaSS
Jxfer = Vxfer * (CaSS - Cai)
    in [mM/ms]
Vxfer = 0.0038 [1/ms]
    in [1/ms]

#
# Calcium dynamics
#
# The equations for the ddt terms were updated in [1]
#
[calcium]
use jup.Jup, jrel.Jrel, jxfer.Jxfer
use phys.F, cell.Cm, cell.Vc, cell.Vss, cell.Vsr
# Free calcium in cytosol, SS, and SR
dot(Cai) = ddt_Cai_total * f_JCai_free
    in [mM]
dot(CaSS) = ddt_CaSS_total * f_JCaSS_free
    in [mM]
dot(CaSR) = ddt_CaSR_total * f_JCaSR_free
    in [mM]
# Derivative of total calcium concentrations (free + buffered)
ddt_Cai_total = -(icab.ICab + ipca.IpCa - 2 * inaca.INaCa) * Cm / (2 * Vc * F) - Jup * Vsr / Vc + Jxfer
    in [mM/ms]
ddt_CaSS_total = -ical.ICaL * Cm / (2 * Vss * F) + Jrel * Vsr / Vss - Jxfer * Vc / Vss
    in [mM/ms]
ddt_CaSR_total = Jup - Jrel
    in [mM/ms]
# Conversion factors from d/dt total to d/dt free
f_JCai_free = 1 / (1 + Buf_c * K_buf_c / (Cai + K_buf_c)^2)
f_JCaSS_free = 1 / (1 + Buf_SS * K_buf_SS / (CaSS + K_buf_SS)^2)
f_JCaSR_free = 1 / (1 + Buf_SR * K_buf_SR / (CaSR + K_buf_SR)^2)
# Buffering: Ca_buffered = Ca_i * Buf_c / (Ca_i + K_bufc)
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]

#
# Sodium dynamics
#
[sodium]
use phys.F, cell.Cm, cell.Vc
dot(Nai) = -INa_total * Cm / (Vc * F)
    in [mM]
INa_total = ina.INa + inab.INab + 3 * inak.INaK + 3 * inaca.INaCa
    in [A/F]

#
# Potassium dynamics
#
[potassium]
use phys.F, cell.Cm, cell.Vc
dot(Ki) = -IK_total * Cm / (Vc * F)
    in [mM]
IK_total = ik1.IK1 + ito.Ito + ikr.IKr + iks.IKs + ipk.IpK + ikur.IKur + stimulus.i_stim - 2 * inak.INaK
    in [A/F]

[[protocol]]
# Level  Start    Length   Period   Multiplier
1        50       0.5      1000.0   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
d = s.run(1000)

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