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

The Na + /K + ATPase is an essential component of cardiac electrophysiology, maintaining physiological Na + and K + concentrations over successive heart beats. Terkildsen et al. (2007) developed a model of the ventricular myocyte Na + /K + ATPase to study extracellular potassium accumulation during ischaemia, demonstrating the ability to recapitulate a wide range of experimental data, but unfortunately there was no archived code associated with the original manuscript. Here we detail an updated version of the model and provide CellML and MATLAB code to ensure reproducibility and reusability. We note some errors within the original formulation which have been corrected to ensure that the model is thermodynamically consistent, and although this required some reparameterisation, the resulting model still provides a good ﬁt to experimental measurements that demonstrate the dependence of Na + /K + ATPase pumping rate upon membrane voltage and metabolite concentrations. To demonstrate thermodynamic consistency we also developed a bond graph version of the model. We hope that these models will be useful for community eﬀorts to assemble a whole-cell cardiomyocyte model which facilitates the investigation of cellular energetics.


Introduction
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., 2007;Terkildsen, 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 and 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 and Gadsby, 1989), extracellular sodium (Nakao and Gadsby, 1989), intracellular sodium (Hansen et al., 2002), extracellular potassium (Nakao and 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., 1971;Gawthrop and Crampin, 2014), and we refer readers to Gawthrop and Smith (1996); Borutzky (2010); and Gawthrop and 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.

Modifications
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 and 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 and 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 ∆G 0 MgATP is the standard free energy of MgATP hydrolysis at pH 0. It appears  Figure 1: Reaction scheme of the cardiac Na + /K + ATPase model. Numbers for each pump state (blue boxes) and reaction names (green) are labelled, with corrected parameters shown in red.
that Terkildsen et al. (2007) started with a standard free energy of −29.6kJ/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 ∆G 0 MgATP = −30.2kJ/mol, which is inconsistent with the typical standard free energy of 11.9kJ/mol at 311K (Tran et al., 2009;Guynn and 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 ∆G 0 MgATP = 11.9kJ/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 and 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 (α + 1 , α + 3 , α − 2 and α − 4 ) were incorrect, leading to inaccurate representations of pump kinetics. We have corrected the equations for these modified rate constants: and ∆ is the unit of charge translocated by the final sodium binding reaction R5. We derive an expression for α + 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".

Reparameterisation of the model
Using the updated model's equations, we fitted parameters to data from Terkildsen et al. (Terkildsen et al., 2007;Terkildsen, 2006). After setting ∆G 0 MgATP to 11.9kJ/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 and Gadsby (1989) was increased from 6× to 15× to obtain a reasonable fit at physiological concentrations.
2. To ensure that cycling velocity magnitudes matched Nakao and Gadsby (1989), the curve for [Na] e = 150mM 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 and 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 and Gadsby, 1989). Updated model parameters are given in Table 1 of Appendix A.
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 and Rudy, 2000;Luo and 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 and Rudy, 2000;Luo and 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  (Faber and Rudy, 2000); (b) The Na + /K + ATPase currents of the original and updated models; (c) A comparison of scaled versions of the updated and original models against the Na + /K + ATPase model in Faber and Rudy (2000). The pump density was increased by a factor of 3.4 in the updated model, and by a factor of 3.9 in the original model.

Bond 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 and Crampin, 2014). The structure of the bond graph model is given in Figure 6 of Appendix B. 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 κ and species thermodynamic constant K (Gawthrop and 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 Ln is the element-wise logarithm operator. The vector k contains the kinetic parameters, λ contains the bond graph parameters, and M contains stoichiometric information. The partitions of these matrices are defined as: where k + = column vector of forward rate constants (12) k − = column vector of reverse rate constants (13) κ = column vector of reaction rate constants (14) K = column vector of species thermodynamic constants (15) 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 (10) can be solved, one possible solution is given by where Exp is the element-wise exponential operator and M † is the Moore-Penrose pseudo-inverse of 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 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 = 38pL was used for the species Na + i , K + i , MgATP, MgADP, Pi and H + , an extracellular volume of W e = 5.182pL was used for Na + e , 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 [Na + ] e = 150mM, 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.

Conclusion
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.