- Author:
- overleaf <overleaf@localhost>
- Date:
- 2017-11-02 00:40:27+00:00
- Desc:
- Uploaded edited paper
- Permanent Source URI:
- https://models.physiomeproject.org/workspace/578/rawfile/e37aa16cf464b16fd113fa6e9612fa5052692b76/Terkildsen_NaK_pump_update.tex
\RequirePackage[l2tabu, orthodox]{nag}
\documentclass[11pt]{article}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{enumitem}
\usepackage[framed,numbered,autolinebreaks,useliterate]{mcode}
\usepackage{fancyvrb}
\usepackage[a4paper,top=3cm]{geometry}
\usepackage{color}
\usepackage{titlesec}
\usepackage{booktabs}
\usepackage{float}
\usepackage{todonotes}
\usepackage{setspace}
\usepackage{fixltx2e}
\usepackage{siunitx}
\usepackage{empheq}
\usepackage[labelfont=bf]{caption}
\usepackage{epstopdf}
\usepackage{empheq}
\usepackage{natbib}
\usepackage{natbibspacing}
\usepackage{etoolbox}
\usepackage[utf8]{inputenc}
\usepackage{datetime}
\usepackage{multirow}
% Use Latin Modern Font
\usepackage{lmodern}
\usepackage[T1]{fontenc}
% Microtype
\usepackage{microtype}
\usepackage[pdftex,bookmarks]{hyperref}
\usepackage{boxedminipage}
%\setlength{\fboxrule}{1pt}
%\setlength{\fboxsep}{12pt}
% Set colours for hyperlinks
\hypersetup{colorlinks=true, linkcolor=blue, citecolor=blue, urlcolor=blue}
% Adjust line spacing for lists
\setlist{itemsep=0pt,leftmargin=*}
% Define a command for reducing size of subscripts
\newcommand{\sss}[1]{\scriptscriptstyle{#1}}
% Change section labelling
\renewcommand{\sectionautorefname}{Section}
\renewcommand{\subsectionautorefname}{{\S}}
\renewcommand{\subsubsectionautorefname}{{\S}}
% Title spacing commands
\titlespacing{\section}{0pt}{0pt}{3pt}
\titlespacing{\subsection}{0pt}{0pt}{0pt}
\setlength{\parindent}{0cm}
\pagestyle{myheadings}
\setlength{\parskip}{0.8\baselineskip}
\newcounter{num}
\setcounter{page}{1}
\titlespacing{\section}{0pt}{*0}{*0}
\newenvironment{customboxedpage}
{\setlength{\fboxrule}{1pt}
\setlength{\fboxsep}{12pt}
\begin{center}
\begin{boxedminipage}{0.7\linewidth}
}{
\end{boxedminipage}
\end{center}
}
% Include comma in citations
\makeatletter
\pretocmd{\NAT@citex}{%
\let\NAT@hyper@\NAT@hyper@citex
\def\NAT@postnote{#2}%
\setcounter{NAT@total@cites}{0}%
\setcounter{NAT@count@cites}{0}%
\forcsvlist{\stepcounter{NAT@total@cites}\@gobble}{#3}}{}{}
\newcounter{NAT@total@cites}
\newcounter{NAT@count@cites}
\def\NAT@postnote{}
% include postnote and \citet closing bracket in hyperlink
\def\NAT@hyper@citex#1{%
\stepcounter{NAT@count@cites}%
\hyper@natlinkstart{\@citeb\@extra@b@citeb}#1%
\ifnumequal{\value{NAT@count@cites}}{\value{NAT@total@cites}}
{\ifNAT@swa\else\if*\NAT@postnote*\else%
\NAT@cmt\NAT@postnote\global\def\NAT@postnote{}\fi\fi}{}%
\ifNAT@swa\else\if\relax\NAT@date\relax
\else\NAT@@close\global\let\NAT@nm\@empty\fi\fi% avoid compact citations
\hyper@natlinkend}
\renewcommand\hyper@natlinkbreak[2]{#1}
% avoid extraneous postnotes, closing brackets
\patchcmd{\NAT@citex}
{\ifNAT@swa\else\if*#2*\else\NAT@cmt#2\fi
\if\relax\NAT@date\relax\else\NAT@@close\fi\fi}{}{}{}
\patchcmd{\NAT@citex}
{\if\relax\NAT@date\relax\NAT@def@citea\else\NAT@def@citea@close\fi}
{\if\relax\NAT@date\relax\NAT@def@citea\else\NAT@def@citea@space\fi}{}{}
\makeatother
\markboth{}{}
\begin{document}
\begin{center}
{\LARGE\textbf{A model of the cardiac Na\textsuperscript{+}/K\textsuperscript{+} ATPase: An update of the Terkildsen model}} \\
Michael Pan, Peter J. Gawthrop, Joseph Cursons, Kenneth Tran, Edmund J. Crampin \\
\longdate\today
\end{center}
\textbf{Abstract} \\
\citet{terkildsen_balance_2007} developed a model of the Na$^+$/K$^+$ ATPase in ventricular myocytes to study extracellular potassium accumulation in ischaemia. The model accounts for a wide range of experimental data, but there is no publicly available code to reproduce the figures in their paper. Furthermore, due to errors in their analysis, the model is not thermodynamically consistent as originally claimed. In this paper we update the model by Terkildsen et al. so that it is reproducible and thermodynamically consistent. The updated model was reparameterised, and provides a good fit to experimental measurements describing the dependence of the pumping rate on membrane voltage and metabolite concentrations. Since the updated model is thermodynamically consistent, it is well-suited for incorporation into a whole-cell model of a cardiomyocyte for the purposes of studying energetics. We developed a bond graph version of the model to demonstrate its thermodynamic consistency, and to contribute towards a thermodynamically consistent whole-cell model of a cardiomyocyte. CellML and MATLAB code has been included to enhance reproducibility and reusability of our model.
\section{Introduction}
In cardiomyocytes, physiological concentrations of Na$^+$ and K$^+$ are maintained by 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. \citet{terkildsen_balance_2007} published a model of the Na$^+$/K$^+$ ATPase in cardiomyocytes to account for the biophysics and thermodynamics of ion transport, with further details described in \citet{terkildsen_modelling_2006}. This model was incorporated into a whole-cell model of a cardiomyocyte to investigate the mechanisms of extracellular potassium accumulation during ischaemia \citep{terkildsen_balance_2007}. Terkildsen et al. found that reduced Na$^+$/K$^+$ ATPase activity was the dominant mechanism for extracellular potassium accumulation. In this paper, we refer to the Na$^+$/K$^+$ ATPase model described in \citet{terkildsen_balance_2007} as the Terkildsen et al. model.
The Na$^+$/K$^+$ ATPase model presented in \citet{terkildsen_balance_2007} was based the model by \citet{smith_development_2004}, and follows the methods proposed by \citet{smith_development_2004} to incorporate thermodynamic constraints and a lumping scheme to simplify the model. A key advantage of the Terkildsen et al. model over the \citet{smith_development_2004} model is that it was fitted to a wider range of data, which included the dependence of pump current on membrane voltage \citep{nakao_[na]_1989}, extracellular sodium \citep{nakao_[na]_1989}, intracellular sodium \citep{hansen_dependence_2002}, extracellular potassium \citep{nakao_[na]_1989} and MgATP \citep{friedrich_na+k+-atpase_1996}. However, the figures for cycling velocity in the original paper are not reproducible using information supplied in the figure legends \citep[Fig. 2]{terkildsen_balance_2007}, and the code used to generate those figures is not publicly available. These reproducibility issues are exacerbated by the use of incorrect equations and numerical values (further described in \autoref{sec:modifications}) which cause physical and thermodynamic inconsistencies.
In this paper, we aim to address the reproducibility and physical issues of the Terkildsen et al. model by updating the model so that it is reproducible and thermodynamically consistent. We modified some equations of the model to ensure thermodynamic consistency (\autoref{sec:modifications}), and refitted the updated model to the same data (\autoref{sec:raparameterisation}). To verify the physical plausibility of the updated model, we developed a bond graph \citep{oster_network_1971,gawthrop_energy-based_2014} version of the updated model. We refer readers to \citet{gawthrop_metamodelling:_1996}; \citet{borutzky_bond_2010}; and \citet{gawthrop_bond-graph_2007} for further information on bond graph theory. Due to its thermodynamic consistency, our updated model is well-suited for incorporation into a thermodynamic model of a cardiomyocyte to study cardiac energetics. MATLAB and CellML \citep{lloyd_cellml:_2004} code has been provided with this paper for reproducibility.
\section{Modifications}
\label{sec:modifications}
The Terkildsen model uses the Post-Albers cycle \citep{apell_electrogenic_1989}, which is a mechanism by which sodium and potassium ions bind one-by-one on one side of the membrane, and unbind on the other side (\autoref{fig:NaK_scheme}). The full Post-Albers cycle was simplified using the methods of \citet{smith_development_2004} to reduce computational complexity. First, the faster reactions were assumed to be in rapid equilibrium to reduce the full 15-state model to a four-state model with eight modified rate constants. The entire cycle was then assumed to be in steady state to further reduce the model to a single equation for the cycle flux, with metabolite dependence incorporated in a manner that accounted for thermodynamic constraints. We found three issues while reimplementing the Terkildsen et al. model, and made several modifications to remedy these issues as described below:
\textbf{Issue 1:} The equilibrium constants are inconsistent with the number of binding sites. Typically, the kinetic rate constants for identical binding sites are assumed to be proportional to the number of binding sites available for binding or unbinding \citep{keener_mathematical_2009}. Thus we modified the reaction scheme in \citet{terkildsen_balance_2007} to be consistent with this assumption (see red parameters in \autoref{fig:NaK_scheme}).
\begin{figure}
\centering
\includegraphics[width=\linewidth]{Figures/Terkildsen_NaK_network}
\caption{\textbf{Reaction scheme of the Terkildsen et al. model.} The numbers for each pump state are labelled in blue boxes, and reaction names are shown in green. Corrected parameters are coloured in red.}
\label{fig:NaK_scheme}
\end{figure}
\textbf{Issue 2:} \citet{terkildsen_balance_2007} derived a detailed balance constraint which was applied during their fitting procedure. This constraint is an equation that relates the kinetic constants defined in \autoref{fig:NaK_scheme}:
\begin{align}
\frac{k_1^+ k_2^+ k_3^+ k_4^+ K_{d,\text{Na}_e}^0(K_{d,\text{Na}_e})^2 (K_{d,\text{K}_i})^2}
{k_1^- k_2^- k_3^- k_4^- K_{d,\text{Na}_i}^0 (K_{d,\text{Na}_i})^2 (K_{d,\text{K}_e})^2 K_{d,\text{MgATP}}} = \exp\left( -\frac{\Delta G_\text{MgATP}^0}{RT} \right)
\label{eq:detailed_balance}
\end{align}
where $R$ 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. \citet{terkildsen_balance_2007} start with a standard free energy of $-29.6\si{kJ/mol}$ at pH 7, but adjust to a physiological pH rather than pH 0. As a result, substituting the model parameter values into equation \eqref{eq:detailed_balance} 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 \citep{tran_thermodynamic_2009,guynn_equilibrium_1973}. Since this results in an overall equilibrium constant over $10^7$-fold greater than the correct value at a temperature of 310K, we decided to use the correct value of $\Delta G_\mathrm{MgATP}^0 = 11.9\si{kJ/mol}$ in the detailed balance constraint.
\textbf{Issue 3:} \citet{terkildsen_balance_2007} used a lumping scheme to reduce the 15-state model to a 4-state model with modified kinetic constants. \citet{terkildsen_balance_2007} use similar expressions for the modified rate constants to \citet{smith_development_2004}, but those expressions were not applicable to the updated assignment of electrical dependence in \citet{terkildsen_balance_2007}. As a result, the expressions for some of the modified rate constants ($\alpha_1^+$, $\alpha_3^+$, $\alpha_2^-$ and $\alpha_4^-$) were incorrect and the original model is an inaccurate representation of pump kinetics. In our updated model, we corrected the equations for these modified rate constants:
\begin{align}
\alpha_1^+ &= \frac{k_1^+ \tilde{\text{Na}}_{i,1}\tilde{\text{Na}}_{i,2}^2}{\tilde{\text{Na}}_{i,1}\tilde{\text{Na}}_{i,2}^2 + (1 + \tilde{\text{Na}}_{i,2})^2 + (1 + \tilde{\text{K}}_i)^2 -1} \\
\alpha_3^+ &= \frac{k_3^+ \tilde{\text{K}}_e^{2}}{\tilde{\text{Na}}_{e,1}\tilde{\text{Na}}_{e,2}^2 + (1 + \tilde{\text{Na}}_{e,2})^2 + (1 + \tilde{\text{K}}_e)^2 -1 } \\
\alpha_2^-&= \frac{k_2^- \tilde{\text{Na}}_{e,1}\tilde{\text{Na}}_{e,2}^2}{\tilde{\text{Na}}_{e,1}\tilde{\text{Na}}_{e,2}^2 + (1 + \tilde{\text{Na}}_{e,2})^2 + (1 + \tilde{\text{K}}_e)^2 -1} \\
\alpha_4^- &= \frac{k_4^- \tilde{\text{K}}_i^2}{\tilde{\text{Na}}_{i,1}\tilde{\text{Na}}_{i,2}^2 + (1 + \tilde{\text{Na}}_{i,2})^2 + (1 + \tilde{\text{K}}_i)^2 -1}
\end{align}
where
\begin{align}
&\tilde{\text{Na}}_{i,1} = \frac{\mathrm{[Na^+]_i}}{K_{d,\text{Na}_i}^0 e^{\Delta FV/RT}}
&&\tilde{\text{Na}}_{i,2} = \frac{\mathrm{[Na^+]_i}}{K_{d,\text{Na}_i}} \\
&\tilde{\text{Na}}_{e,1} = \frac{\mathrm{[Na^+]_e}}{K_{d,\text{Na}_e^0} e^{(1+\Delta) zFV/RT}}
&&\tilde{\text{Na}}_{e,2} = \frac{\mathrm{[Na^+]_e}}{K_{d,\text{Na}_e}} \\
&\tilde{\text{K}}_i = \frac{\mathrm{[K^+]_i}}{K_{d,\text{K}_i}}
&&\tilde{\text{K}}_e = \frac{\mathrm{[K^+]_e}}{K_{d,\text{K}_e}}
\end{align}
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$,
\begin{align}
\alpha_1^+ &= k_1^+ \frac{x_6}{x_6 + x_5 + x_4 + x_3 + x_2 + x_1} \notag \\
&= k_1^+ \frac{1}{1 + x_5/x_6 + x_4/x_6 + x_3/x_6 + x_2/x_6 + x_1/x_6} \notag \\
&= \frac{k_1^+}{1 + 2\tilde{\text{Na}}_{i,1}^{-1} + 2\tilde{\text{Na}}_{i,1}^{-1}\tilde{\text{Na}}_{i,2}^{-1} + \tilde{\text{Na}}_{i,1}^{-1}\tilde{\text{Na}}_{i,2}^{-2} + 2\tilde{\text{Na}}_{i,1}^{-1}\tilde{\text{Na}}_{i,2}^{-2}\tilde{\text{K}}_i + \tilde{\text{Na}}_{i,1}^{-1}\tilde{\text{Na}}_{i,2}^{-2}\tilde{\text{K}}_i^2} \notag \\
&= \frac{k_1^+ \tilde{\text{Na}}_{i,1}\tilde{\text{Na}}_{i,2}^2}{\tilde{\text{Na}}_{i,1}\tilde{\text{Na}}_{i,2}^2 + (1 + \tilde{\text{Na}}_{i,2})^2 + (1 + \tilde{\text{K}}_i)^2 -1}
\end{align}
Because it was not possible to fix the above issues without significantly changing the kinetics of the model, we reparameterised 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 \citep{terkildsen_balance_2007} as the ``original model''.
\section{Reparameterisation of the model}
\label{sec:raparameterisation}
Using the updated model's equations, we fitted parameters to data from Terkildsen et al. \citep{terkildsen_balance_2007,terkildsen_modelling_2006}. The original model was parameterised by minimising an objective function that measured the divergence of the model from the data. We parameterised the updated model using the methods of \citet{terkildsen_modelling_2006}, and setting $\Delta G_\mathrm{MgATP}^0$ to its correct value of $11.9\si{kJ/mol}$ in the thermodynamic constraint (Equation \eqref{eq:detailed_balance}). Other minor changes to the fitting procedure include:
\begin{enumerate}
\item The weighting for extracellular potassium above 5.4 mM for the data of \citet{nakao_[na]_1989} was increased from $6\times$ to $15\times$ to obtain a reasonable fit at physiological concentrations.
\item To ensure that the resulting cycling velocities had magnitudes that matched \citet{nakao_[na]_1989}, the curve for $[\mathrm{Na}]_e=150 \si{mM}$ was fitted in its un-normalised form.
\item Rather than using a local optimiser with literature sources for initial parameter estimates, we used particle swarm optimisation \citep{kennedy_particle_1995} followed by a local optimiser to find a global minimum of the objective function.
\end{enumerate}
The results of the fitting process are shown in Figures \ref{fig:fitting} and \ref{fig:metabolite_dependence}. The updated model provides good fits to each data source, and the quality of fits are comparable to \citet{terkildsen_modelling_2006}. Despite a slightly worse fit at lower extracellular sodium concentrations compared to the original model (\autoref{fig:fitting}(a)) the unnormalised cycling velocities (\autoref{fig:fitting}(b)) appear to be more consistent with experimental data that suggest saturated cycling velocity at positive membrane potentials is relatively insensitive to extracellular sodium \citep{nakao_[na]_1989}. Updated model parameters are given in \autoref{tab:Terkildsen_parameters} of \autoref{sec:parameters}.
\begin{figure}
\centering
\begin{tabular}{c c}
(a) & (b) \\
\includegraphics[width=0.4\linewidth]{Figures/Terkildsen_kinetic_fit_KA_comparison} &
\includegraphics[width=0.4\linewidth]{Figures/Terkildsen_vcyc} \\
(c) & \multirow{2}{0.4\linewidth}[0.9cm]{
\begin{minipage}{\linewidth}
\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*}
\end{minipage} } \\
\includegraphics[width=0.4\linewidth]{Figures/Terkildsen_fit_NG_I} &
\end{tabular}
\caption{\textbf{Fit of the updated Terkildsen model to current-voltage measurements.} \textbf{(a)} Comparison of model to extracellular sodium and voltage data \citep[Fig. 3]{nakao_[na]_1989}. Cycling velocities were normalised to a value of $55\si{s^{-1}}$ at $V = 40\si{mV}$. \textbf{(b)} Unnormalised cycling velocties. \textbf{(c)} Comparison of model to whole-cell current measurements \citep[Fig. 2(a)]{nakao_[na]_1989}.}
\label{fig:fitting}
\end{figure}
% * <ktra014@aucklanduni.ac.nz> 2017-09-22T02:43:11.081Z:
%
% Figure 2a shows unnormalised cycling velocity (same as 2b) - should it be normalised?
%
% ^ <panm@student.unimelb.edu.au> 2017-09-25T01:00:04.320Z:
%
% The cycling velocity is "normalised", although I'm not sure if that's the right word as it's done in a different way from 2d-f. In Figure 2a, the curves were scaled so that they all passed through 55s^-1 at -40mV - which is how the model was compared to data in the original paper. I plot the cycling velocities without this scaling in Figure 2b, although it turns out there scaling didn't change the values by much. The scaling had a more noticeable impact in the original paper.
%
% ^.
\begin{figure}
\centering
\begin{tabular}{c c}
(a) & \multirow{2}{0.4\linewidth}[1cm]{
\begin{minipage}{\linewidth}
\small
\begin{align*}
V &= -40\si{mV}\\
\mathrm{[Na^+]_e} &= 140\si{mM}\\
\mathrm{[K^+]_i} &= 70\si{mM}\\
\mathrm{[K^+]_e} &= 5.4\si{mM}\\
\mathrm{pH} &= 7.2\\
\mathrm{[Pi]_{tot}} &= 1\si{mM}\\
\mathrm{[MgATP]} &= 2\si{mM}\\
\mathrm{[MgADP]} &= 0\si{mM}\\
T &= 309\si{K}
\end{align*}
\end{minipage}}\\
\includegraphics[width=0.4\linewidth]{Figures/Terkildsen_fit_Hansen_Nai_comparison} & \\[0.2cm]
(b) & \multirow{2}{0.4\linewidth}[1cm]{
\begin{minipage}{\linewidth}
\small
\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*}
\end{minipage}}\\
\includegraphics[width=0.4\linewidth]{Figures/Terkildsen_fit_NG_Ke_comparison} & \\[0.2cm]
(c) & \multirow{2}{0.4\linewidth}[1cm]{
\begin{minipage}{\linewidth}
\small
\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*}
\end{minipage}}\\
\includegraphics[width=0.4\linewidth]{Figures/Terkildsen_fit_Friedrich_MgATP_comparison} &
\end{tabular}
\caption{\textbf{Fit of the updated Terkildsen model to metabolite dependence data.} Simulation conditions are displayed on the right of each figure. \textbf{(a)} Comparison of model to intracellular sodium data \citep[Fig. 7(a)]{hansen_dependence_2002}, normalised to the cycling velocity at $\mathrm{[Na]_i = 50\si{mM}}$. \textbf{(b)} Comparison of model to extracellular potassium data \citep[Fig. 11(a)]{nakao_[na]_1989}, normalised to the cycling velocity at $\mathrm{[K]_e = 5.4\si{mM}}$. \textbf{(c) } Comparison of model to ATP data \citep[Fig. 3(b)]{friedrich_na+k+-atpase_1996}, normalised to the cycling velocity at $\mathrm{[MgATP] = 0.6\si{mM}}$.}
\label{fig:metabolite_dependence}
\end{figure}
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 (LRd00) model \citep{faber_action_2000,luo_dynamic_1994} (\autoref{fig:updated_original_comp}(a)). A comparison of the response between the updated model and the original model is shown in \autoref{fig:updated_original_comp}(b). The two models behave almost identically at resting membrane potentials, but the updated model has a much higher current during the action potential. As noted in \citet{terkildsen_modelling_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 LRd00 model. We compared scaled versions of the updated and original models to the Na$^+$/K$^+$ ATPase current described using equations in the LRd00 model and found that both versions of the Terkildsen et al. model exhibit similar behaviour to the description in the LRd00 model (\autoref{fig:updated_original_comp}(c)). The updated model behaves more closely to the LRd00 Na$^+$/K$^+$ ATPase model because it has a more variable current. Thus we hypothesise that under physiological concentrations, the updated model has better compatibility with the LRd00 whole-cell model. A CellML version of the updated model is included with this paper, and reproduces the curve for $\mathrm{[Na]_e} = 150\si{mM}$ in \autoref{fig:fitting}(a).
% * <ktra014@aucklanduni.ac.nz> 2017-09-22T03:20:41.581Z:
%
% I wonder if its worth using SED-ML to create settings for recreating all the figures. But I assume your matlab code does all that already?
%
% ^ <panm@student.unimelb.edu.au> 2017-09-25T01:30:52.320Z:
%
% My MATLAB code can recreate the figures (although I haven't tested it on other computers). I haven't learned how to write SED-ML yet, but will look into it.
%
% ^.
\begin{figure}
\centering
\begin{tabular}{c c c}
(a) & (b) & (c) \\
\includegraphics[width=0.3\linewidth]{Figures/LRd_AP_waveform} &
\includegraphics[width=0.3\linewidth]{Figures/NaK_current} &
\includegraphics[width=0.3\linewidth]{Figures/NaK_current_scaled}
\end{tabular}
\caption{\textbf{A comparison of the updated Terkildsen model to existing models.} (a) The action potential waveform used to simulate the pumps \citep{faber_action_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 \citet{faber_action_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. $\mathrm{[Na^+]_i} = 10\si{mM}$, $\mathrm{[Na^+]_e} = 140\si{mM}$, $\mathrm{[K^+]_i} = 145\si{mM}$, $\mathrm{[K^+]_e} = 5.4\si{mM}$, $\mathrm{pH} = 7.095$, $\mathrm{[Pi]_{tot}} = 0.8\si{mM}$, $\mathrm{[MgATP]} = 6.95\si{mM}$, $\mathrm{[MgADP]} = 0.035\si{mM}$, $T = 310\si{K}$.}
\label{fig:updated_original_comp}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=0.7\linewidth]{Figures/CellML_comp}
\caption{\textbf{A comparison of the kinetic and bond graph models.} $\mathrm{[Na^+]_i} = 50\si{mM}$, $\mathrm{[Na^+]_e} = 150\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}$. In the bond graph model, zero concentrations were approximated by using a concentration of 0.001mM to avoid undefined values.}
\label{fig:kinetic_BG_comp}
\end{figure}
\section{Bond graph model}
To verify the physical plausibility of the updated model, and to aid its incorporation into larger models, we developed a bond graph version of the updated model. The bond graph model represents the full unsimplified biochemical cycle, and reactions that were assumed to be in rapid equilibrium were replaced by reactions with fast kinetic parameters with the same equilibrium constant. The structure of the bond graph is given in \autoref{fig:Terkildsen_NaK} of \autoref{sec:BG_structure}. We find the parameters of the bond graph model \citep{gawthrop_hierarchical_2015} by solving the matrix equation
\begin{align}
\textbf{Ln}(\mathbf{k}) = \mathbf{M} \textbf{Ln}(\mathbf{W} \boldsymbol{\lambda}) \label{eq:bg_general}
\end{align}
where $\textbf{Ln}$ is the element-wise logarithm operator and
\begin{align}
\textbf{k} = \begin{bmatrix}
k^+ \\ k^- \\ k^c
\end{bmatrix}, \quad
\textbf{M} = \left[ \begin{array}{c | c}
I_{n_r \times n_r} & {N^f}^T \\ \hline
I_{n_r \times n_r} & {N^r}^T \\ \hline
0 & N^c
\end{array} \right], \quad
\boldsymbol{\lambda} = \begin{bmatrix}
\kappa \\ K
\end{bmatrix}
\label{eq:standard}
\end{align}
We define the partitions of $\textbf{k}$, $\boldsymbol{\lambda}$ and $\mathbf{M}$ as
\begin{align}
k^+ &= \text{column vector of forward rate constants} \\
k^- &= \text{column vector of reverse rate constants} \\
\kappa &= \text{column vector of reaction rate constants} \\
K &= \text{column vector of species thermodynamic constants} \\
N^f &= \text{forward stoichiometric matrix} \\
N^r &= \text{reverse stoichiometric matrix}
\end{align}
The matrices $k^c$ and $N^c$ were used to enforce further constraints between species thermodynamic constants. These constraints describe the equilibria of the identical ions in different compartments, and the equilibrium constant of ATP hydrolysis. Assuming that equation \eqref{eq:bg_general} can be solved, one possible solution is given by
\begin{align}
\boldsymbol{\lambda_0 } = \mathbf{W}^{-1} \textbf{Exp} (\mathbf{M}^\dagger \textbf{Ln} (\mathbf{k}))
\end{align}
where $\textbf{Exp}$ is the element-wise exponential operator and $\mathbf{M}^\dagger$ is the Moore-Penrose pseudoinverse of $\mathbf{M}$. For this bond graph model, we use the diagonal matrix $\mathbf{W}$ to scale the species thermodynamic constants according to the volume they reside in. For consistency with \citet{terkildsen_balance_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 \autoref{fig:fitting}(a) to reproduce the curve for $\mathrm{[Na]_e}=150\si{mM}$. The membrane voltage was increased at a slow rate to induce quasi-steady-state behaviour. There is a high degree of correspondence between the kinetic and bond graph models (\autoref{fig:kinetic_BG_comp}). The only substantial difference occurs near the initial voltage of $V=-120\si{mV}$, where the bond graph model moves towards its steady state. The CellML code describing the bond graph model is provided with this paper, and reproduces the curve in \autoref{fig:kinetic_BG_comp}.
\section{Conclusion}
In this paper, we updated the Na$^+$/K$^+$ ATPase model by \citet{terkildsen_balance_2007} to make it reproducible and thermodynamically consistent. We found errors in the analysis in the original paper which resulted in inaccurate descriptions of pump kinetics and ultimately a thermodynamically inconsistent model. Thus the equations and constraints within the model were corrected to make it thermodynamically consistent. We refitted the updated model to data, and simulations matched experimental measurements well. Since the updated model is thermodynamically consistent, it has a natural bond graph representation. We developed a bond graph version of the model to demonstrate thermodynamic consistency. CellML and MATLAB code for both the kinetic and bond graph models have been included with this paper to aid reproducibility. The thermodynamic consistency and reusability of our updated model make it ideal for incorporation into future whole-cell models to study cardiac cell energetics.
\bibliographystyle{model2-names_edit}
\small
\bibliography{bibliography}{}
\normalsize
\appendix
\section{Parameters}
\label{sec:parameters}
\begin{table}[H]
\caption{\textbf{Kinetic parameters for the updated Terkildsen et al. model.} Refer to \autoref{fig:NaK_scheme} for a schematic.}
\centering
\bgroup
\def\arraystretch{1.3}
\begin{tabular}{c p{0.48\linewidth} l}
\toprule
\textbf{Parameter} & \textbf{Description}& \textbf{Value} \\ \midrule
$k_1^+$ & Forward rate constant of reaction R6 & $1423.2\ \si{s^{-1}}$\\
$k_1^-$ & Reverse rate constant of reaction R6 & $225.9048\ \si{s^{-1}}$\\
$k_2^+$ & Forward rate constant of reaction R7 & $11564.8064\ \si{s^{-1}}$\\
$k_2^-$ & Reverse rate constant of reaction R7 & $36355.3201\ \si{s^{-1}}$\\
$k_3^+$ & Forward rate constant of reaction R13 & $194.4506\ \si{s^{-1}}$\\
$k_3^-$ & Reverse rate constant of reaction R13 & $281037.2758\ \si{mM^{-2}s^{-1}}$\\
$k_4^+$ & Forward rate constant of reaction R15 & $30629.8836\ \si{s^{-1}}$\\
$k_4^-$ & Reverse rate constant of reaction R15 & $1.574\times 10^{6}\ \si{s^{-1}}$\\
$K_{\text{d,Nai}}^0$ & Voltage-dependent dissociation constant of intracellular $\mathrm{Na^+}$ & $579.7295\ \si{mM}$\\
$K_{\text{d,Nae}}^0$ & Voltage-dependent dissociation constant of extracellular $\mathrm{Na^+}$ & $0.034879\ \si{mM}$\\
$K_{\text{d,Nai}}$ & Voltage-independent dissociation constant of intracellular $\mathrm{Na^+}$ & $5.6399\ \si{mM}$\\
$K_{\text{d,Nae}}$ & Voltage-independent dissociation constant of extracellular $\mathrm{Na^+}$ & $10616.9377\ \si{mM}$\\
$K_{\text{d,Ki}}$ & Dissociation constant of intracellular $\mathrm{K^+}$ & $16794.976\ \si{mM}$\\
$K_{\text{d,Ke}}$ & Dissociation constant of extracellular $\mathrm{K^+}$ & $1.0817\ \si{mM}$\\
$K_{\text{d,MgATP}}$ & Dissociation constant of MgATP & $140.3709\ \si{mM}$\\
$\Delta$ & Charge translocated by reaction R5 & -0.0550 \\
Pump density & Number of pumps per $\si{\mu m^2}$ & $1360.2624\ \si{\mu m^{-2}}$ \\ \bottomrule&
\end{tabular}
\egroup
\label{tab:Terkildsen_parameters}
\end{table}
\begin{table}[H]
\caption{\textbf{Parameters for the bond graph version of the updated Terkildsen et al. model.} Parameters were derived by using an intracellular volume of 38pL and an extracellular volume of 5.182pL. Refer to \autoref{fig:Terkildsen_NaK} for the bond graph schematic.}
\centering
\begin{tabular}{cl c l}
\toprule
\textbf{Component} & \textbf{Description}& \textbf{Parameter} & \textbf{Value} \\ \midrule
R1 & Reaction R1 & $\kappa_1$ & $ 330.5462\ \si{fmol/s} $\\
R2 & Reaction R2 & $\kappa_2$ & $ 132850.9145\ \si{fmol/s} $\\
R3 & Reaction R3 & $\kappa_3$ & $ 200356.0223\ \si{fmol/s} $\\
R4 & Reaction R4 & $\kappa_4$ & $ 2238785.3951\ \si{fmol/s} $\\
R5 & Reaction R5 & $\kappa_5$ & $ 10787.9052\ \si{fmol/s} $\\
R6 & Reaction R6 & $\kappa_6$ & $ 15.3533\ \si{fmol/s} $\\
R7 & Reaction R7 & $\kappa_7$ & $ 2.3822\ \si{fmol/s} $\\
R8 & Reaction R8 & $\kappa_8$ & $ 2.2855\ \si{fmol/s} $\\
R9 & Reaction R9 & $\kappa_9$ & $ 1540.1349\ \si{fmol/s} $\\
R10 & Reaction R10 & $\kappa_{10}$ & $ 259461.6507\ \si{fmol/s} $\\
R11 & Reaction R11 & $\kappa_{11}$ & $ 172042.3334\ \si{fmol/s} $\\
R12 & Reaction R12 & $\kappa_{12}$ & $ 6646440.3909\ \si{fmol/s} $\\
R13 & Reaction R13 & $\kappa_{13}$ & $ 597.4136\ \si{fmol/s} $\\
R14 & Reaction R14 & $\kappa_{14}$ & $ 70.9823\ \si{fmol/s} $\\
R15 & Reaction R15 & $\kappa_{15}$ & $ 0.015489\ \si{fmol/s} $\\
$\text{P}_1$ & Pump state ATP--E\textsubscript{i} K\textsubscript{2}
& $K_1$ & $101619537.2009\ \si{fmol^{-1}}$ \\
$\text{P}_2$ & Pump state ATP--E\textsubscript{i} K\textsubscript{1}
& $K_2$ & $63209.8623\ \si{fmol^{-1}}$ \\
$\text{P}_3$ & Pump state ATP--E\textsubscript{i}
& $K_3$ & $157.2724\ \si{fmol^{-1}}$ \\
$\text{P}_4$ & Pump state ATP--E\textsubscript{i} Na\textsubscript{1}
& $K_4$ & $14.0748\ \si{fmol^{-1}}$ \\
$\text{P}_5$ & Pump state ATP--E\textsubscript{i} Na\textsubscript{2}
& $K_5$ & $5.0384\ \si{fmol^{-1}}$ \\
$\text{P}_6$ & Pump state ATP--E\textsubscript{i} Na\textsubscript{3}
& $K_6$ & $92.6964\ \si{fmol^{-1}}$ \\
$\text{P}_7$ & Pump state P--E\textsubscript{i} (Na\textsubscript{3})
& $K_7$ & $4854.5924\ \si{fmol^{-1}}$ \\
$\text{P}_8$ & Pump state P--E\textsubscript{e} Na\textsubscript{3}
& $K_8$ & $15260.9786\ \si{fmol^{-1}}$ \\
$\text{P}_9$ & Pump state P--E\textsubscript{e} Na\textsubscript{2}
& $K_9$ & $13787022.8009\ \si{fmol^{-1}}$ \\
$\text{P}_{10}$ & Pump state P--E\textsubscript{e} Na\textsubscript{1}
& $K_{10}$ & $20459.5509\ \si{fmol^{-1}}$ \\
$\text{P}_{11}$ & Pump state P--E\textsubscript{e}
& $K_{11}$ & $121.4456\ \si{fmol^{-1}}$ \\
$\text{P}_{12}$ & Pump state P--E\textsubscript{e} K\textsubscript{1}
& $K_{12}$ & $3.1436\ \si{fmol^{-1}}$ \\
$\text{P}_{13}$ & Pump state P--E\textsubscript{e} K\textsubscript{2}
& $K_{13}$ & $0.32549\ \si{fmol^{-1}}$ \\
$\text{P}_{14}$ & Pump state E\textsubscript{e} (K\textsubscript{2})
& $K_{14}$ & $156.3283\ \si{fmol^{-1}}$ \\
$\text{P}_{15}$ & Pump state ATP--Ee (K\textsubscript{2})
& $K_{15}$ & $1977546.8577\ \si{fmol^{-1}}$ \\
Ki & Intracellular $\text{K}_\text{i}^+$ & $K_\text{Ki}$ & $0.0012595\ \si{fmol^{-1}}$ \\
Ke & Extracellular $\text{K}_\text{e}^+$ & $K_\text{Ke}$ & $0.009236\ \si{fmol^{-1}}$ \\
Nai & Intracellular $\text{Na}_\text{i}^+$ & $K_\text{Nai}$ & $0.00083514\ \si{fmol^{-1}}$ \\
Nae & Extracellular $\text{Na}_\text{e}^+$ & $K_\text{Nae}$ & $0.0061242\ \si{fmol^{-1}}$ \\
$\text{MgATP}$ & Intracellular MgATP & $K_\text{MgATP}$ & $2.3715\ \si{fmol^{-1}}$ \\
$\text{MgADP}$ & Intracellular MgADP & $K_\text{MgADP}$ & $7.976 \times 10^{-5} \ \si{fmol^{-1}}$ \\
$\text{P}_\text{i}$ & Free inorganic phosphate & $K_\text{Pi}$ & $0.04565\ \si{fmol^{-1}}$ \\
H & Intracellular $\text{H}^+$ & $K_\text{H}$ & $0.04565\ \si{fmol^{-1}}$ \\
mem & Membrane capacitance & $C_m$ & $153400\ \si{fF}$ \\
zF{\_}5 & Charge translocated by R5 & $z_5$ & $-0.0550$ \\
zF{\_}8 & Charge translocated by R8 & $z_8$ & $-0.9450$ \\ \bottomrule& &
\end{tabular}
\label{tab:bg_parameters}
\end{table}
\newpage
\section{Bond graph model structure}
\label{sec:BG_structure}
\begin{figure}[H]
\centering
\includegraphics[width=\linewidth]{Figures/Terkildsen_NaK_BG_full}
\caption{\textbf{Bond graph structure of the Terkildsen et al. model.} Pump states are coloured in blue, and reactions are coloured in green. The names for these components are consistent with their labels in \autoref{fig:NaK_scheme}.}
\label{fig:Terkildsen_NaK}
\end{figure}
\end{document}
