The cardiac Na⁺/K⁺ ATPase: An updated, thermodynamically consistent model

1Introduction¶

Cardiomyocytes maintain Na$^+$ and K$^+$ ions within their physiological concentration range, in part by using Na$^+$/K$^+$ ATPases located on their plasma membranes. The Na$^+$/K$^+$ ATPase is an electrogenic ion pump that uses energy from ATP hydrolysis to drive the transport of Na$^+$ and K$^+$ ions against an electrochemical gradient. A previous model of the cardiomyocyte Na$^+$/K$^+$ ATPase was developed by Terkildsen et al. (2007) and subsequently incorporated into a whole-cell cardiomyocyte model to demonstrate that reduced Na$^+$/K$^+$ ATPase activity plays a dominant role in extracellular potassium accumulation during ischaemia Terkildsen et al., 2007Terkildsen, 2006. In this manuscript the Na$^+$/K$^+$ ATPase model of Terkildsen et al. (2007) will subsequently be referred to as the Terkildsen et al. model.

The Na$^+$/K$^+$ ATPase model presented in Terkildsen et al. (2007) was based upon an earlier implementation which proposed thermodynamic constraints and a lumping scheme for model simplification Smith & Crampin, 2004. A key development was that the Terkildsen et al. model could reproduce a wider range of data which captured the dependence of the pump current upon membrane voltage Nakao & Gadsby, 1989, extracellular sodium Nakao & Gadsby, 1989, intracellular sodium Hansen et al., 2002, extracellular potassium Nakao & Gadsby, 1989 and MgATP Friedrich et al., 1996. Unfortunately the cycling velocity figures presented within the original paper are not reproducible using information supplied in the figure legends Terkildsen et al., 2007, Fig. 2 and code used to generate those figures was not publicly archived. These issues are exacerbated by apparent errors within the reported equations and parameter values (further described in Section 2) which result in physical and thermodynamic inconsistencies.

Here we address these issues, updating the model to ensure that it is thermodynamically consistent, and archiving MATLAB and CellML Lloyd et al., 2004 code for reproducibility. This required the modification of several equations (Section 2) and re-parameterisation through fitting to the original data (Section 3). To verify the physical plausibility of the updated model we have also developed a bond graph version Oster et al., 1971Gawthrop & Crampin, 2014, and we refer readers to Gawthrop & Smith (1996); Borutzky (2010); and Gawthrop & Bevan (2007) for further information on bond graph theory. Given the thermodynamic consistency of our updated model we believe that it is particularly well-suited for incorporation into community efforts for developing a thermodynamic model of a cardiomyocyte to ultimately study whole-heart cardiac energetics.

2Modifications¶

The Terkildsen et al. model uses the Post-Albers cycle Apell, 1989, a model in which sodium and potassium ions bind individually on one side of the membrane, and unbind on the other side (Figure 1). The full Post-Albers cycle was simplified to reduce computational complexity by assuming that faster reactions are in rapid equilibrium, reducing the full 15-state model to a four-state model with eight modified rate constants Smith & Crampin, 2004. The entire cycle was then assumed to be in steady-state such that the model simplified to a single equation for cycling velocity, with metabolite dependence incorporated in a manner that accounted for thermodynamic constraints. We identified three issues while reimplementing the Terkildsen et al. model, and made several modifications to remedy these issues:

Issue 1: Equilibrium constants were inconsistent with the number of binding sites. For identical binding sites the kinetic rate constants are typically assumed to be proportional to the number of sites available for binding/unbinding Keener & Sneyd, 2009 and we modified the reaction scheme from Terkildsen et al. (2007) to achieve this (red parameters within Figure 1).

Issue 2: The detailed balance constraint used during fitting procedure appears to have used an incorrect parameter value with important consequences on the thermodynamic consistency of the model. This constraint relates the kinetic constants defined in Figure 1:

where $R=8.314$ J/mol/K is the universal gas constant, $T$ is the absolute temperature, and $\Delta G_\mathrm{MgATP}^0$ is the standard free energy of MgATP hydrolysis at pH 0. It appears that Terkildsen et al. (2007) started with a standard free energy of $-29.6\si{kJ/mol}$ at pH 7, but adjusted to a physiological pH rather than pH 0. As a result, substituting the model parameter values into equation (1) results in $\Delta G_\mathrm{MgATP}^0 = -30.2\si{kJ/mol}$, which is inconsistent with the typical standard free energy of $11.9\si{kJ/mol}$ at 311K Tran et al., 2009Guynn & Veech, 1973. At a temperature of 310K this results in an overall equilibrium constant over $10^7$-fold greater than the correct value. Thus we use $\Delta G_\mathrm{MgATP}^0 = 11.9\si{kJ/mol}$ within the detailed balance constraint.

Issue 3: The lumping scheme used to reduce the 15-state model to a 4-state model with modified kinetic constants was similar to Smith & Crampin (2004), however with the updated assignment of electrical dependence in Terkildsen et al. (2007) some expressions were not applicable. As a result, expressions for several modified rate constants ($\alpha_1^+$, $\alpha_3^+$, $\alpha_2^-$ and $\alpha_4^-$) were incorrect, leading to inaccurate representations of pump kinetics. We have corrected the equations for these modified rate constants:

where

and $\Delta$ is the unit of charge translocated by the final sodium binding reaction R5. We derive an expression for $\alpha_1^+$, and expressions for the other modified rate constants follow similarly. Since the pump states 1 to 6 are lumped together, the constant $k_1^+$ is scaled according to the ratio between the amount of state 6 and the total amount of states 1–6. If we represent $x_i$ as the molar amount of state $i$, then:

Because it was not possible to fix the above issues without significantly changing the kinetics of the model, we subsequently had to reparameterise the Terkildsen et al. model such that it would be physically and thermodynamically consistent. In subsequent sections, we shall refer to the reparameterised model with updated equations as the “updated model” and the model with equations and parameters described in Terkildsen et al. (2007) as the “original model”.

3Reparameterisation of the model¶

Using the updated model’s equations, we fitted parameters to data from Terkildsen et al. (2007)Terkildsen (2006). After setting $\Delta G_\mathrm{MgATP}^0$ to $11.9\si{kJ/mol}$ in the thermodynamic constraint, Equation (1), we parameterised the updated model using similar methods to the original model Terkildsen, 2006 which minimised an objective function describing divergence of the model from experimental data. Minor changes to the fitting procedure include:

1. The weighting for extracellular potassium above 5.4 mM for the data of Nakao & Gadsby (1989) was increased from $6\times$ to $15\times$ to obtain a reasonable fit at physiological concentrations.

2. To ensure that cycling velocity magnitudes matched Nakao & Gadsby (1989), the curve for $[\mathrm{Na}]_e=150 \si{mM}$ was fitted without normalisation.

3. Rather than using a local optimiser with literature sources for initial parameter estimates, we minimised the objective function by using particle swarm optimisation Kennedy & Eberhart, 1995 followed by a local optimiser to find a global minimum.

The updated model provides good fits to each data source (Figures 2 & 3), and the quality of fits are comparable to Terkildsen (2006), although we achieved a slightly worse fit at lower extracellular sodium concentrations (Figure 2(A)). It should be noted, however, that our model appears to be more consistent with experimental data that suggest saturated cycling velocity at positive membrane potentials is relatively insensitive to extracellular sodium (Figure 2(B)) Nakao & Gadsby, 1989. Updated model parameters are given in Table 1 of Section 8.

A	$$\begin{align*} \mathrm{[Na^+]_i} &= 50\si{mM} \\ \mathrm{[K^+]_i} &= 0\si{mM} \\ \mathrm{[K^+]_e} &= 5.4\si{mM} \\ \mathrm{pH} &= 7.4 \\ \mathrm{[Pi]_{tot}} &= 0\si{mM} \\ \mathrm{[MgATP]} &= 10\si{mM} \\ \mathrm{[MgADP]} &= 0\si{mM} \\ T &= 310\si{K} \end{align*}$$ (5)
B

A	$$\begin{align*} V &= 0\si{mV}\\ \mathrm{[Na^+]_e} &= 0\si{mM}\\ \mathrm{[K^+]_i} &= 80\si{mM}\\ \mathrm{[K^+]_e} &= 15\si{mM}\\ \mathrm{pH} &= 7.2\\ \mathrm{[Pi]_{tot}} &= 1\si{mM}\\ \mathrm{[MgATP]} &= 2\si{mM}\\ \mathrm{[MgADP]} &= 0\si{mM}\\ T &= 308\si{K} \end{align*}$$ (6)
B	$$\begin{align*} V &= 0\si{mV}\\ \mathrm{[Na^+]_i} &= 50\si{mM}\\ \mathrm{[Na^+]_e} &= 150\si{mM}\\ \mathrm{[K^+]_i} &= 140\si{mM}\\ \mathrm{pH} &= 7.4\\ \mathrm{[Pi]_{tot}} &= 0.5\si{mM}\\ \mathrm{[MgATP]} &= 10\si{mM}\\ \mathrm{[MgADP]} &= 0.02\si{mM}\\ T &= 310\si{K} \end{align*}$$ (7)
C	$$\begin{align*} V &= 0\si{mV}\\ \mathrm{[Na^+]_i} &= 40\si{mM}\\ \mathrm{[Na^+]_e} &= 0\si{mM}\\ \mathrm{[K^+]_i} &= 0\si{mM}\\ \mathrm{[K^+]_e} &= 5\si{mM}\\ \mathrm{pH} &= 7.4\\ \mathrm{[Pi]_{tot}} &= 0\si{mM}\\ \mathrm{[MgADP]} &= 0\si{mM}\\ T &= 297\si{K} \end{align*}$$ (8)

The response of the updated model to an action potential input was simulated by using an action potential waveform generated from the Luo and Rudy (2000) model Faber & Rudy, 2000Luo & Rudy, 1994 (Figure 4(A)). The original and updated models behave almost identically at resting membrane potentials, but the updated model has a much higher current during the action potential (Figure 4(B)). As noted in Terkildsen (2006), the current of the pump is far lower at physiological intracellular sodium concentrations, thus the pump density needs to be appropriately scaled to be compatible with the Luo-Rudy model. Scaled versions of the Na$^+$/K$^+$ ATPase current within the updated and original models are qualitatively similar to that described using the Luo-Rudy equations Faber & Rudy, 2000Luo & Rudy, 1994, however there are some differences in the resulting waveforms (Figure 4(C)). In particular, the updated model behaves more similarly to the Luo-Rudy Na$^+$/K$^+$ ATPase formulation because it has a more variable current, and thus we hypothesise that under physiological concentrations, the updated model is more compatible with the whole-cell model by Luo and Rudy. A CellML version of the updated model is included with this manuscript.

A	B	C

4Bond graph model¶

To verify the physical plausibility of the updated model and to aid incorporation into larger models of cardiac energetics, we have also developed a bond graph version. Bond graphs are an energy-based approach to modelling physical systems, thus they ensure thermodynamic consistency Gawthrop & Crampin, 2014. The structure of the bond graph model is given in Figure 6 of Section 9. The process of converting the model into a bond graph required two notable changes to its representation. Firstly, the bond graph model represents the full unsimplified biochemical cycle, and reactions originally assumed to be in rapid equilibrium were replaced by reactions with fast kinetic parameters that conferred the same equilibrium constant. Thus the bond graph model contains 15 states, and is a close but not exact approximation of the kinetic model. Secondly, because kinetic parameters are often thermodynamically inconsistent Liebermeister et al., 2010, the bond graph approach requires chemical reaction networks to be specified using a different set of parameters: the reaction rate constant $\kappa$ and species thermodynamic constant $K$ Gawthrop & Crampin, 2014. These parameters always describe thermodynamically consistent systems, regardless of their numerical value. As a result, we converted the kinetic parameters of our model into an equivalent set of bond graph parameters Gawthrop et al., 2015 by using the following matrix equation:

where $\textbf{Ln}$ is the element-wise logarithm operator. The vector $\mathbf{k}$ contains the kinetic parameters, $\boldsymbol{\lambda}$ contains the bond graph parameters, and $\mathbf{M}$ contains stoichiometric information. The partitions of these matrices are defined as:

where

The matrices $k^c$ and $N^c$ were used to enforce further constraints between the thermodynamic constants of different species, in particular the equilibria of individual ions present within different compartments, and the equilibrium constant of ATP hydrolysis. Assuming that equation (9) can be solved, one possible solution is given by

where $\textbf{Exp}$ is the element-wise exponential operator and $\mathbf{M}^\dagger$ is the Moore-Penrose pseudo-inverse of $\mathbf{M}$. Since the bond graph framework is energy-based, species must be expressed as molar amounts rather than concentrations to adequately compare energies in different compartments. Therefore we use the diagonal matrix $\mathbf{W}$ to scale the species thermodynamic constants according to the volume of the compartments they reside in. For consistency with Terkildsen et al. (2007), an intracellular volume of $W_i = 38\si{pL}$ was used for the species $\mathrm{Na_i^+}$, $\mathrm{K_i^+}$, $\mathrm{MgATP}$, $\mathrm{MgADP}$, $\mathrm{Pi}$ and $\mathrm{H^+}$, an extracellular volume of $W_e = 5.182\si{pL}$ was used for $\mathrm{Na_e^+}$, $\mathrm{K_e^+}$, and a constant of 1 was used for each of the pump states.

The bond graph model was simulated under the conditions described in Figure 2(A) to reproduce the curve for $\mathrm{[Na^+]_e}=150\si{mM}$, using a slowly increasing membrane voltage to induce quasi-steady-state behaviour. There is a high degree of correspondence between the kinetic and bond graph models (Figure 5). The closeness of the two models suggests that the fast kinetic constants are good approximations of reactions in rapid equilibrium, although we note that the deviation between the bond graph and kinetic models increases slightly at higher cycling velocities when the faster reactions begin to limit flux through the cycle. CellML code describing the bond graph model and reproducing the curve in Figure 5 is provided with this manuscript.

5Conclusion¶

In this manuscript we describe an updated model for the cardiac Na$^+$/K$^+$ ATPase model originally developed by Terkildsen et al. (2007). We have corrected errors with the original model formulation and refitted necessary parameters to ensure that the resulting model is thermodynamically consistent while still recapitulating a wide range of experimental data. We note that the updated model has a natural bond graph representation, and include CellML and MATLAB code for both the kinetic and bond graph models to aid reproducibility. We believe that the thermodynamic consistency and improved reusability of our updated model make it ideal for incorporation into future whole-cell models to study cardiac cell energetics.

6Instructions for running code¶

6.1MATLAB¶

The MATLAB code is located in the MATLAB/code directory. Figures 2–3 can be generated by running the script fit_NaK_parameters.m. Figure 4 can be generated by running NaK_sim.m, and Figure 5 can be generated by running CellML_model_comparison.m.

6.2CellML¶

The figures within this paper can be generated by opening the SED-ML files in OpenCOR. The corresponding figures are referenced within the file names.

8Parameters¶

9Bond graph model structure¶