From 9ea4811b43588b359f49d2442f8ab1e6ab0875f7 Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Fri, 16 Aug 2024 16:11:13 +0100 Subject: [PATCH] Added Fabbri 2017 human SAN model --- c/fabbri-2017.mmt | 627 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 627 insertions(+) create mode 100644 c/fabbri-2017.mmt diff --git a/c/fabbri-2017.mmt b/c/fabbri-2017.mmt new file mode 100644 index 0000000..d2679f6 --- /dev/null +++ b/c/fabbri-2017.mmt @@ -0,0 +1,627 @@ +[[model]] +author: Michael Clerx, Aditi Agrawal +name: fabbri-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] +C = 5.7e-5 [uF] + in [uF] +dot(V) = -i_tot / C + in [mV] +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] + +# +# 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] +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) / (k * (1 + iso_slope / 100)))) + k = 4.3371 [mV] + in [mV] + iso_shift = piecewise(iso.iso > 0, -8 [mV], 0 [mV]) + in [mV] + iso_slope = piecewise(iso.iso > 0, -27, 0) + tau = 0.001 / (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/s/mV] * (av + 41.8 [mV]) / (exp(-(av + 41.8 [mV]) / 2.5 [mV]) - 1) - 0.0849 [1/mV/s] * (av + 6.8 [mV]) / (exp(-(av + 6.8 [mV]) / 4.8 [mV]) - 1) + in [1/s] + beta = 0.01143 [1/s/mV] * (bv + 1.8 [mV]) / (exp((bv + 1.8 [mV]) / 2.5 [mV]) - 1) + in [1/s] +dot(f) = (inf - f) / tau + tau = 0.001 * (44.3 [s] + 230 [s] * 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 * inf / alpha + in [s] + Km = 0.000338 [mM] + in [mM] + alpha = 0.0075 [1/s] + in [1/s] +# 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 * 0.6 * (65.17 / (0.57 [1/s] * exp(-0.08 [1/mV] * (V + 44 [mV])) + 0.065 [1/s] * exp(0.1 [1/mV] * (V + 45.93 [mV]))) + 10.1 [s]) + in [s] +dot(r) = (inf - r) / tau + inf = 1 / (1 + exp(-(V - 19.3 [mV]) / 15 [mV])) + tau = 0.001 * 0.66 * 1.4 * (15.59 [s] / (1.037 * exp(0.09 [1/mV] * (V + 30.61 [mV])) + 0.369 * exp(-0.12 [1/mV] * (V + 23.84 [mV]))) + 2.98 [s]) + in [s] + +# +# Rapid delayed-rectifier potassium current +# +[ikr] +use membrane.V +gKr = 0.00424 [uS] + in [uS] +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 extra.Na_o, extra.Ca_o +use membrane.V +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 + calcium.Ca_sub / Kci * (1 + exp(-Qci * V / phys.RTF) + sodium.Na_i / Kcni) + sodium.Na_i / K1ni * (1 + sodium.Na_i / K2ni * (1 + sodium.Na_i / K3ni)) +do = 1 + Ca_o / Kco * (1 + exp(Qco * V / phys.RTF)) + extra.Na_o / K1no * (1 + Na_o / K2no * (1 + Na_o / K3no)) +k12 = calcium.Ca_sub / Kci * exp(-Qci * V / phys.RTF) / di +k21 = Ca_o / Kco * exp(Qco * V / phys.RTF) / do +k23 = Na_o / K1no * Na_o / K2no * (1 + Na_o / K3no) * exp(-Qn * V / (2 * phys.RTF)) / do +k32 = exp(Qn * V / (2 * phys.RTF)) +k34 = Na_o / (K3no + Na_o) +k43 = sodium.Na_i / (K3ni + sodium.Na_i) +k14 = sodium.Na_i / K1ni * sodium.Na_i / K2ni * (1 + sodium.Na_i / K3ni) * exp(Qn * V / (2 * phys.RTF)) / di +k41 = exp(-Qn * V / (2 * phys.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] +use cell.vi, cell.vsub +dot(Na_i) = (1 - clamp) * -INa_tot / ((vi + 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() +