INaK fixes in O'Hara model

c/ohara-2011.mmt

@@ -7,13 +7,36 @@ desc: """
     The 2011 model for the undiseased human ventricular AP and calcium
     transient by O'Hara et al. [1].
 
-    The model builds on several other models, including [2].
+    The model builds on several other models, including [2]. 
 
     This implementation is based on the original Matlab code provided by the
     authors. It was checked against the original code by comparing the
-    calculated derivatives. Where possible, the page of the appendix to [1]
-    that describes each component is indicated.
+    calculated derivatives. After this, two small changes to the INaK 
+    formulation were made (see below).
 
+    Comments and errata on [1]:
+
+    1. The comment "Publisher's Note: Errors in Text S1" describes three issues
+       in the supplement text, which are not present in the published model
+       code. Similarly, "Publisher's Note: Error in Figure 1, panel C, 5th 
+       subplot" and "Publisher's Note: Rate dependent APD curves for TP model
+       were arranged incorrectly in Figure 18A" describe issues in the figures.
+       No changes to the Myokit model were made based on these comments.
+    2. The comment "A Note About the Fast Na Current (I~~Na~~) Formulation"
+       explains issues with using the O'Hara model for tissue simulations, and
+       suggests a replacement of the INa formulation for this purpose. No
+       changes to the Myokit model were made based on this comment.
+    4. The comment "Ko and liquid junction potential" describes what could be a
+       model issue (no LJP correction), but gives no solution, and so no
+       changes to the Myokit model were made based on this note.
+    5. The comment "Two minor issues in the INaK model implementation"
+       describes two errors in the INaK equations, which do not significantly
+       affect model output under baseline conditions. Both are corrected in
+       this Myokit model.
+
+    Where possible, the page of the supplement to [1] that describes each 
+    component has been indicated.
+
     References:
 
     [1] O'Hara, T., Virág, L., Varró, A., & Rudy, Y. (2011). Simulation of the
@@ -741,8 +764,11 @@ INaCa_ss = 0.2 * inaca.fNaCa * inaca.gNaCa * allo * (JncxNa + 2 * JncxCa)
 
 #
 # INaK: Sodium/potassium ATPase current
-# Based on Smith and Crampin, 2004, PBMB
-# Page 14
+# 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
+#
+# Supplement page 14
 #
 [inak]
 use membrane.V
@@ -779,8 +805,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]
@@ -811,22 +838,21 @@ 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]