Added ORd-cipa INaK corrections

myokit · Mar 7, 2024 · b894d0f · b894d0f
1 parent 28d3738
commit b894d0f
Showing 1 changed file with 59 additions and 52 deletions.
diff --git a/c/ohara-cipa-v1-2017.mmt b/c/ohara-cipa-v1-2017.mmt
@@ -1,6 +1,6 @@
 [[model]]
 name: li-2017
-version: 20230208
+version: 20240307
 mmt_authors: Michael Clerx
 display_name: O'Hara-CiPA-v1, 2017
 desc: """
@@ -18,6 +18,10 @@ desc: """
      - Conductances for INaL, ICaL, IKs, and IK1 were rescaled, as described
        in [6].
 
+    This implementation includes a further two corrections to INaK, as 
+    documented at
+    https://docs.google.com/document/d/111fqNzQGvGAjB_PrkvejEhzqwROrR6czz_OBz7Ep-iM
+
     References:
 
     [1] Li, Z., Dutta, S., Sheng, J., Tran, P. N., Wu, W., Chang, K., Mdluli,
@@ -114,7 +118,7 @@ pace = 0 bind pace
 
 #
 # Membrane potential
-# Page 5 of the appendix to [2]
+# Page 5 of the supplement to [2]
 #
 [membrane]
 dot(V) = -(i_ion + stimulus.i_stim)
@@ -133,7 +137,7 @@ i_ion = (
 
 #
 # Stimulus current
-# Page 5 of the appendix to [2]
+# Page 5 of the supplement to [2]
 #
 [stimulus]
 i_stim = engine.pace * amplitude
@@ -143,11 +147,11 @@ amplitude = -80 [A/F]
 
 #
 # Cell geometry
-# Page 6 of the appendix to [2]
+# Page 6 of the supplement to [2]
 #
 [cell]
 mode = 1
-    desc: The type of cell. Endo = 0, Epi = 1, Mid = 2
+    desc: The type of cell. Endo = 0, Epi = 1, M = 2
 L = 0.01 [cm] : Cell length
     in [cm]
 r = 0.0011 [cm] : Cell radius
@@ -196,7 +200,7 @@ FFRT = F * FRT
 
 #
 # Extracellular concentrations
-# Page 5 of the appendix to [2]
+# Page 5 of the supplement to [2]
 #
 [extra]
 Na_o = 140 [mM] : Extracellular Na+ concentration
@@ -208,7 +212,7 @@ K_o = 5.4 [mM] : Extracellular K+ concentration
 
 #
 # Reversal potentials
-# Page 6 of the appendix to [2]
+# Page 6 of the supplement to [2]
 #
 [rev]
 ENa = phys.RTF * log(extra.Na_o / sodium.Na_i)
@@ -220,12 +224,12 @@ EK = phys.RTF * log(extra.K_o / potassium.K_i)
 PNaK = 0.01833
     desc: Permeability ratio K+ to Na+
 EKs = phys.RTF * log((extra.K_o + PNaK * extra.Na_o) / (potassium.K_i + PNaK * sodium.Na_i))
-    desc: Reversal potential for IKs
     in [mV]
+    desc: Reversal potential for IKs
 
 #
 # INa: Fast sodium current
-# Page 6 of the appendix to [2]
+# Page 6 of the supplement to [2]
 #
 # The fast sodium current is modelled using a Hodgkin-Huxley type formulation
 # including activation (m), slow and fast components of inactivation (h) and
@@ -239,19 +243,19 @@ use membrane.V
 sm = 1 / (1 + exp((V + 39.57 [mV]) / -9.871 [mV]))
     desc: Steady state value for m-gates
 tm = 1 [ms] / (6.765 * exp((V + 11.64 [mV]) / 34.77 [mV]) + 8.552 * exp(-(V + 77.42 [mV]) / 5.955 [mV]))
-    desc: Time constant for m-gates
     in [ms]
+    desc: Time constant for m-gates
 dot(m) = (sm - m) / tm
     desc: Activation gate for INa
 # h-gates
 sh = 1 / (1 + exp((V + 82.9 [mV]) / 6.086 [mV]))
     desc: Steady-state value for h-gates
 thf = 1 [ms] / (1.432e-5 * exp((V + 1.196 [mV]) / -6.285 [mV]) + 6.149 * exp((V + 0.5096 [mV]) / 20.27 [mV]))
-    desc: Time constant for fast development of inactivation in INa
     in [ms]
+    desc: Time constant for fast development of inactivation in INa
 ths = 1 [ms] / (0.009794 * exp((V + 17.95 [mV]) / -28.05 [mV]) + 0.3343 * exp((V + 5.73 [mV]) / 56.66 [mV]))
-    desc: Time constant for slow development of inactivation in INa
     in [ms]
+    desc: Time constant for slow development of inactivation in INa
 Ahf = 0.99 : Fraction of INa channels with fast inactivation
 Ahs = 1 - Ahf : Fraction of INa channels with slow inactivation
 dot(hf) = (sh - hf) / thf
@@ -264,14 +268,14 @@ h = Ahf * hf + Ahs * hs
 sj = sh
     desc: Steady-state value for j-gate in INa
 tj = 2.038 [ms] + 1 [ms] / (0.02136 * exp((V + 100.6 [mV]) / -8.281 [mV]) + 0.3052 * exp((V + 0.9941 [mV]) / 38.45 [mV]))
-    desc: Time constant for j-gate in INa
     in [ms]
+    desc: Time constant for j-gate in INa
 dot(j) = (sj - j) / tj
     desc: Recovery from inactivation gate for non-phosphorylated INa
 # Phosphorylated channels
 thsp = 3 * ths
-    desc: Time constant for h-gate of phosphorylated INa
     in [ms]
+    desc: Time constant for h-gate of phosphorylated INa
 shsp = 1 / (1 + exp((V + 89.1 [mV]) / 6.086 [mV]))
     desc: Steady-state value for h-gate of phosphorylated INa
 dot(hsp) = (shsp - hsp) / thsp
@@ -292,7 +296,7 @@ INa = gNa * (V - rev.ENa) * m^3 * ((1 - camk.f) * h * j + camk.f * hp * jp)
 
 #
 # INaL: Late component of the sodium current
-# Page 7 of the appendix to [2]
+# Page 7 of the supplement to [2]
 #
 [inal]
 use membrane.V
@@ -324,15 +328,15 @@ INaL = fNaL * gNaL * (V - rev.ENa) * m * ((1 - camk.f) * h + camk.f * hp)
 
 #
 # Ito: Transient outward potassium current
-# page 8 of the appendix to [2]
+# page 8 of the supplement to [2]
 #
 [ito]
 use membrane.V
 ta = 1.0515 [ms] / (one + two)
+    in [ms]
     one = 1 / (1.2089 * (1 + exp((V - 18.4099 [mV]) / -29.3814 [mV])))
     two = 3.5 / (1 + exp((V + 100 [mV]) / 29.3814 [mV]))
     desc: Time constant for Ito activation
-    in [ms]
 sa = 1 / (1 + exp((V - 14.34 [mV]) / -14.82 [mV]))
     desc: Steady-state value for Ito activation
 dot(a) = (sa - a) / ta
@@ -388,7 +392,7 @@ Ito = fto * gto * (V - rev.EK) * ((1 - camk.f) * a * i + camk.f * ap * ip)
 # ICaL: L-type calcium current
 # ICaNa: Sodium current through the L-type calcium channel
 # ICaK: Potassium current through the L-type calcium channel
-# Page 9 of the appendix to [2]
+# Page 9 of the supplement to [2]
 #
 # The ICaL channel is modeled using activation, inactivation (fast and slow),
 # Ca-dependent inactivation (fast and slow) and recovery from Ca-dependent
@@ -621,7 +625,7 @@ IKr = fKr * gKr * sqrt(extra.K_o / 5.4 [mM]) * O * (V - rev.EK)
 
 #
 # IKs: Slow delayed rectifier potassium current
-# Page 11 of the appendix to [2]
+# Page 11 of the supplement to [2]
 #
 # Modelled with two activation gates
 #
@@ -649,7 +653,7 @@ IKs = fKs * gKs * KsCa * x1 * x2 * (V - rev.EKs)
 
 #
 # IK1: Inward rectifier potassium current
-# Page 12 of the appendix to [2]
+# Page 12 of the supplement to [2]
 #
 # Modelled with an activation channel and an instantaneous inactivation channel
 #
@@ -670,12 +674,12 @@ gK1 = 0.1908 [mS/uF]
     in [mS/uF]
     desc: Maximum conductance for IK1 (before scaling)
 IK1 = fK1 * gK1 * sqrt(K_o / 1 [mM]) * r * x * (V - rev.EK)
-    desc: Inward rectifier Potassium current
     in [A/F]
+    desc: Inward rectifier Potassium current
 
 #
 # INaCa: Sodium/calcium exchange current
-# page 12 of the appendix to [2]
+# Page 12 of the supplement to [2]
 #
 [inaca]
 use membrane.V
@@ -770,7 +774,7 @@ INaCa_total = INaCa + inacass.INaCa_ss
 
 #
 # INaCa_ss: Sodium/calcium exchanger current into the L-type subspace
-# Page 12 of the appendix to [2]
+# Page 12 of the supplement to [2]
 #
 [inacass]
 use membrane.V
@@ -833,13 +837,16 @@ JncxNa = 3 * (E4 * k7 - E1 * k8) + E3 * k4pp - E2 * k3pp
 JncxCa = E2 * k2 - E1 * k1
     in [1/s]
 INaCa_ss = 0.2 * inaca.fNaCa * inaca.gNaCa * allo * (JncxNa + 2 * JncxCa)
-    desc: Sodium/Calcium exchange current into the T-Tubule subspace
     in [A/F]
+    desc: Sodium/Calcium exchange current into the T-Tubule subspace
 
 #
 # INaK: Sodium/potassium ATPase current
-# Based on Smith and Crampin, 2004, PBMB
-# Page 14 of the appendix to [2]
+# Based on Smith & Crampin 2004 https://doi.org/10.1016/j.pbiomolbio.2004.01.010
+# The formulation below was corrected from the O'Hara implementation, see
+# https://docs.google.com/document/d/111fqNzQGvGAjB_PrkvejEhzqwROrR6czz_OBz7Ep-iM
+#
+# Page 14 of the supplement to [2]
 #
 [inak]
 use membrane.V
@@ -876,8 +883,9 @@ MgATP = 9.8 [mM]
     in [mM]
 Kmgatp = 1.698e-7 [mM]
     in [mM]
-H = 1e-7 [mM]
+H = 1e-4 [mM]
     in [mM]
+    note: Corrected to 1e-7 [M] (pH 7) from original value of 1e-7 [mM]
 eP = 4.2 [mM]
     in [mM]
 Khp = 1.698e-7 [mM]
@@ -908,32 +916,31 @@ a4 = k4p * MgATP / Kmgatp / (1 + MgATP / Kmgatp)
     in [Hz]
 b4 = k4m * (K_i / Kki)^2 / ((1 + Na_i / Knai)^3 + (1 + K_i / Kki)^2 - 1)
     in [Hz]
-x1 = a4 * a1 * a2 + b2 * b4 * b3 + a2 * b4 * b3 + b3 * a1 * a2
+x1 = a4 * a1 * a2 + b1 * b4 * b3 + a2 * b4 * b3 + b3 * a1 * a2
     in [Hz^3]
+    note: Corrected from the original code (b1 in second term)
 x2 = b2 * b1 * b4 + a1 * a2 * a3 + a3 * b1 * b4 + a2 * a3 * b4
     in [Hz^3]
 x3 = a2 * a3 * a4 + b3 * b2 * b1 + b2 * b1 * a4 + a3 * a4 * b1
     in [Hz^3]
 x4 = b4 * b3 * b2 + a3 * a4 * a1 + b2 * a4 * a1 + b3 * b2 * a1
     in [Hz^3]
-E1 = x1 / (x1 + x2 + x3 + x4)
-E2 = x2 / (x1 + x2 + x3 + x4)
-E3 = x3 / (x1 + x2 + x3 + x4)
-E4 = x4 / (x1 + x2 + x3 + x4)
-JnakNa = 3 * (E1 * a3 - E2 * b3)
+r = (a1 * a2 * a3 * a4 - b1 * b2 * b3 * b4) / (x1 + x2 + x3 + x4)
     in [1/s]
-JnakK = 2 * (E4 * b1 - E3 * a1)
+JnakNa = 3 * r
     in [1/s]
+JnakK = -2 * r
+    in [1/s]    
 fNaK = piecewise(cell.mode == 0, 1, cell.mode == 1, 0.9, 0.7)
 PNaK = 30 [C/F]
     in [C/F]
 INaK = fNaK * PNaK * (JnakNa + JnakK)
-    desc: Sodium/Potassium ATPase current
     in [A/F]
+    desc: Sodium/Potassium ATPase current
 
 #
 # IKb: Background potassium current
-# Page 15 of the appendix to [2]
+# Page 15 of the supplement to [2]
 #
 [ikb]
 use membrane.V
@@ -942,12 +949,12 @@ fKb = if(cell.mode == 1, 0.6, 1)
 gKb = 0.003 [mS/uF]
     in [mS/uF]
 IKb = fKb * gKb * xkb * (V - rev.EK)
-    desc: Background Potassium current
     in [A/F]
+    desc: Background Potassium current
 
 #
 # INab: Background sodium current
-# Page 15 of the appendix to [2]
+# Page 15 of the supplement to [2]
 #
 [inab]
 use membrane.V
@@ -960,7 +967,7 @@ INab = PNab * V * phys.FFRT * (sodium.Na_i * evf - extra.Na_o) / (evf - 1)
 
 #
 # ICab: Background calcium current
-# Page 15 of the appendix to [2]
+# Page 15 of the supplement to [2]
 #
 [icab]
 use membrane.V
@@ -973,7 +980,7 @@ ICab = PCab * 4 * V * phys.FFRT * (calcium.Ca_i * evf2 - 0.341 * extra.Ca_o) / (
 
 #
 # IpCa: Sarcolemmal calcium pump current
-# Page 15 of the appendix to [2]
+# Page 15 of the supplement to [2]
 #
 [ipca]
 GpCa = 0.0005 [A/F]
@@ -984,7 +991,7 @@ IpCa = GpCa * calcium.Ca_i / (0.0005 [mM] + calcium.Ca_i)
 
 #
 # Jrel: SR Calcium release flux via ryanodine receptor
-# Page 17 of the appendix to [2]
+# Page 17 of the supplement to [2]
 #
 [ryr]
 use calcium.Ca_jsr
@@ -1000,8 +1007,8 @@ dot(Jrelnp) = (inf - Jrelnp) / tau
             in [ms]
     inf = base * if(cell.mode == 2, 1.7, 1)
         in [mM/ms]
-    base = -1 [mM/ms/mV] * a_rel * ical.ICaL / (1 + (1.5 [mM] / Ca_jsr)^8)
-        in [mM/ms]
+        base = -1 [mM/ms/mV] * a_rel * ical.ICaL / (1 + (1.5 [mM] / Ca_jsr)^8)
+            in [mM/ms]
 btp = 1.25 * bt
     in [ms]
 a_relp = 0.5 * btp
@@ -1014,15 +1021,15 @@ dot(Jrelp) = (inf - Jrelp) / tau
             in [ms]
     inf = base * if(cell.mode == 2, 1.7, 1)
         in [mM/ms]
-    base = -1 [mM/ms/mV] * a_relp * ical.ICaL / (1 + (1.5 [mM] / Ca_jsr)^8)
-        in [mM/ms]
+        base = -1 [mM/ms/mV] * a_relp * ical.ICaL / (1 + (1.5 [mM] / Ca_jsr)^8)
+            in [mM/ms]
 Jrel = (1 - camk.f) * Jrelnp + camk.f * Jrelp
     desc: SR Calcium release flux via Ryanodine receptor
     in [mM/ms]
 
 #
 # Jup: Calcium uptake via SERCA pump
-# Page 17 of the appendix to [2]
+# Page 17 of the supplement to [2]
 #
 [serca]
 use calcium.Ca_i, calcium.Ca_jsr, calcium.Ca_nsr
@@ -1042,7 +1049,7 @@ Jtr = (Ca_nsr - Ca_jsr) / 100 [ms]
 
 #
 # Diffusion fluxes
-# Page 16 of the appendix to [2]
+# Page 16 of the supplement to [2]
 #
 [diff]
 JdiffNa = (sodium.Na_ss - sodium.Na_i) / 2 [ms]
@@ -1054,7 +1061,7 @@ Jdiff = (calcium.Ca_ss - calcium.Ca_i) / 0.2 [ms]
 
 #
 # Intracellular sodium concentrations
-# Page 18 of the appendix to [2]
+# Page 18 of the supplement to [2]
 #
 [sodium]
 use cell.AFC, cell.vss, cell.vmyo
@@ -1070,7 +1077,7 @@ dot(Na_ss) = -INa_ss_tot * AFC / vss - diff.JdiffNa
 
 #
 # Intracellular potassium concentrations
-# Page 18 of the appendix to [2]
+# Page 18 of the supplement to [2]
 #
 [potassium]
 use cell.AFC, cell.vss, cell.vmyo
@@ -1094,7 +1101,7 @@ dot(K_ss) = -IK_ss_tot * AFC / vss - diff.JdiffK
 
 #
 # Intracellular calcium concentrations and buffers
-# Page 18 of the appendix to [2]
+# Page 18 of the supplement to [2]
 #
 [calcium]
 use cell.AFC, cell.vmyo, cell.vnsr, cell.vjsr, cell.vss
@@ -1124,13 +1131,13 @@ ICa_tot = ipca.IpCa + icab.ICab - 2 * inaca.INaCa
     in [A/F]
 dot(Ca_i) = buff * (-ICa_tot * AFC / (2 * vmyo) - serca.Jup * vnsr / vmyo + diff.Jdiff * vss / vmyo)
     in [mM]
-    desc: Intracellular Calcium concentratium
+    desc: Intracellular calcium concentration
     buff = 1 / (1 + cmdnmax * kmcmdn / (kmcmdn + Ca_i)^2 + trpnmax * kmtrpn / (kmtrpn + Ca_i)^2)
 ICa_ss_tot = ical.ICaL - 2 * inacass.INaCa_ss
     in [A/F]
 dot(Ca_ss) = buff * (-ICa_ss_tot * AFC / (2 * vss) + ryr.Jrel * vjsr / vss - diff.Jdiff)
     in [mM]
-    desc: Calcium concentratium in the T-Tubule subspace
+    desc: Calcium concentration in the T-Tubule subspace
     buff = 1 / (1 + BSRmax * KmBSR / (KmBSR + Ca_ss)^2 + BSLmax * KmBSL / (KmBSL + Ca_ss)^2)
 dot(Ca_jsr) = buff * (serca.Jtr - ryr.Jrel)
     in [mM]