Location: The cardiac Na+ /K+ ATPase: An updated, thermodynamically consistent model @ 5669cde2344c / Terkildsen_NaK_pump_update.tex

Author:
overleaf <overleaf@localhost>
Date:
2017-11-02 00:40:19+00:00
Desc:
Initial upload
Permanent Source URI:
https://models.physiomeproject.org/workspace/578/rawfile/5669cde2344c8b922964b7c646b4379de7b1149c/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}

% 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}}

% 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{Updating the Terkildsen Na$^+$/K$^+$ ATPase model}} \\
	Michael Pan, Peter J. Gawthrop, Joseph Cursons, Kenneth Tran, Edmund J. Crampin \\
	\longdate\today
\end{center}
\section{Introduction}
Terkildsen et al. developed a model of the Na$^+$/K$^+$ ATPase in cardiomyocytes to account for biophysics and thermodynamics of ion transport \citep{terkildsen_balance_2007,terkildsen_modelling_2006}. The Na$^+$/K$^+$ ATPase is an electrogenic ion pump which uses energy from ATP hydrolysis to power the transport of ions against an electrochemical gradient. The Terkildsen et al. model is based on the \citet{smith_development_2004} Na$^+$/K$^+$ ATPase model, and fits the model to a wider range of data. The Terkildsen et al. model follows the methods proposed by Smith and Crampin to incorporate thermodynamic constraints and simplify the model by using a lumping scheme. The Terkildsen et al. model of the Na$^+$/K$^+$ ATPase 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_modelling_2006}. It was found that reduced Na$^+$/K$^+$ ATPase activity was the dominant mechanism for extracellular potassium accumulation. In this paper, we refer to Na$^+$/K$^+$ ATPase model described in \citet{terkildsen_balance_2007} rather than the whole-cell model as the Terkildsen et al. model.

A key advantage of the Terkildsen et al. model over the earlier model by Smith and Crampin is that it was fitted to a wider range of data. The model was fitted to data from \citet{nakao_[na]_1989}, where the steady-state pump cycling rate was found for a range of membrane voltages and extracellular sodium concentrations. Additional data for whole-cell currents was used to parametrise the pump density of the model. The Na$^+$/K$^+$ ATPase is also activated by intracellular sodium, extracellular potassium, and MgATP. Data describing the activation of the Na$^+$/K$^+$ ATPase with these metabolites was also used to parameterise the model \citep{hansen_dependence_2002,nakao_[na]_1989,friedrich_na+k+-atpase_1996}.

A limitation of the presentation of the model in the original paper is that figures for the cycling velocity of the Na$^+$/K$^+$ ATPase are not reproducible using information from the figure legends \citep[Fig. 2]{terkildsen_balance_2007}, and the code to generate those figures is not publicly available. These reproducibility issues are exacerbated by the fact that there are physical issues in the model that arise from the use of incorrect equations and numerical values (further described in \autoref{sec:modifications}). In this paper, we modify the original model to ensure physical and thermodynamic consistency (\autoref{sec:modifications}), and refit the updated model to data (\autoref{sec:modifications}). To verify the physical plausibility of the updated model, we develop a bond graph \citep{gawthrop_energy-based_2014,borutzky_bond_2010} version of the updated model. We hope that this bond graph model will aid in the incorporation of the model in to a thermodynamic model of a cardiomyocyte. MATLAB and CellML \citep{lloyd_cellml:_2004} code has been provided with this paper to ensure reproducibility of the model and figures.

\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 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, which reduced the original 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 flux, with metabolite dependence incorporated in a manner that accounted for thermodynamic constraints.

We found three issues in the Terkildsen model, and made modifications to fix 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. 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.} Corrected parameters are coloured in red.}
\label{fig:NaK_scheme}
\end{figure}

\textbf{Issue 2:} \citet{terkildsen_balance_2007} apply a detailed balance constraint during their fitting process so that the cycle has zero flux when the metabolites are in equilibrium:
\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}
However, solving for the standard free energy of MgATP hydrolysis $\Delta G_\mathrm{MgATP}$ by substituting the parameter values resulted in a value of $-30.2\si{kJ/mol}$, which is inconsistent with the typical standard free energy of $11.9\si{kJ/mol}$ \citep{tran_thermodynamic_2009,guynn_equilibrium_1973}. Since this results in an overall equilibrium constant over $10^7$-fold greater than the correct value, we decided modify the thermodynamic constraint to use the correct value of $\Delta G_\mathrm{MgATP}$.

\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. However, the equations for some of the modified rate constants ($\alpha_1^+$, $\alpha_3^+$, $\alpha_2^-$ and $\alpha_4^-$) were incorrect. While the incorrect equations do not appear to cause thermodynamic inconsistencies, they are nonetheless an inaccurate representation of pump kinetics. In our updated model, we use the corrected 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}

Because it was not possible to fix the above issues without significantly changing the kinetics of the model, we decided to update and reparametersise 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''.

\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) & (d) \\
		\includegraphics[width=0.4\linewidth]{Figures/Terkildsen_fit_NG_I} &
		\includegraphics[width=0.4\linewidth]{Figures/Terkildsen_fit_Hansen_Nai_comparison} \\
		(e) & (f) \\
		\includegraphics[width=0.4\linewidth]{Figures/Terkildsen_fit_NG_Ke_comparison} &
		\includegraphics[width=0.4\linewidth]{Figures/Terkildsen_fit_Friedrich_MgATP_comparison} 
	\end{tabular}
	\caption{\textbf{Fit of the updated Terkildsen model to data.} \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}$. $\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}$. \textbf{(b)} Unnormalised cycling velocties under the same conditions as (a). \textbf{(c)} Comparison of model to whole-cell current measurements \citep[Fig. 2(a)]{nakao_[na]_1989} under the same conditions as (a). \textbf{(d)} 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}}$. $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}$. \textbf{(e)} 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}}$. $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}$. \textbf{(f) } 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}}$. $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}$.} 
	\label{fig:fitting}
\end{figure}

\section{Reparameterisation of the model}
\label{sec:raparameterisation}
Using the equations of the updated model, we fitted the parameters of the model to the data as in Terkildsen et al. \citep{terkildsen_balance_2007,terkildsen_modelling_2006}. Terkildsen et al. parameterised their model by minimising an objective function that measured the divergence of the model from the data. We parameterised the updated model using the methods of Terkildsen et al. \citep{terkildsen_modelling_2006}, setting $\Delta G_\mathrm{MgATP}$ to its correct value of $11.9\si{kJ/mol}$ in the thermodynamic constraint (Equation \eqref{eq:detailed_balance}). There were some other minor changes to the fitting process:
\begin{enumerate}
	\item The weighting for extracellular potassium above 5.4mM was 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 unnormalised form.
	\item Rather than using literature sources for initial parameter estimates, we used particle swarm optimisation followed by a local optimiser to minimise the objective function.
\end{enumerate}
The results of the fitting process are shown in \autoref{fig:fitting}. The updated model provides good fits to each of the data sources used. The parameters of the updated model are given in \autoref{sec:parameters}.

The updated model was simulated by inputting an action potential waveform generated from the Luo and Rudy 2000 (LRd00) model \citet{faber_action_2000,luo_dynamic_1994} (\autoref{fig:updated_original_comp}(a)). A comparison of the updated model to the original model is shown in \autoref{fig:updated_original_comp}(b). The two models behave very similarly 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. Thus we compared scaled versions of the updated and original Terkildsen et al. models to the Na$^+$/K$^+$ ATPase current described using equations in the LRd00 model. We find that both versions of the Terkildsen et al. model exhibit similar behaviour to the description in the LRd00 model. The updated model has a far greater variability in current across an action potential, and behaves similar to the LRd00 Na$^+$/K$^+$ ATPase. 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).

\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 constants with the same equilibrium constant. The structure of the bond graph is given in Figures \ref{fig:Terkildsen_NaK} and \ref{fig:BG_submodules} 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}
Here $k^+$ and $k^-$ are column vectors of the forward and reverse rate constants respectively, at zero membrane voltage. $\kappa$ and $K$ are column vectors of the reaction rate constants and species thermodynamic constants respectively. $N^f$ and $N^r$ are the forward and reverse stoichiometric matrices respectively. The matrices $k^c$ and $N^c$ were used to enforce further constraints between species thermodynamic constants. These constraints describe the similarity of the same ions in different compartments, and the equilibrium constant of ATP hydrolysis. Assuming that equation \eqref{eq:bg_general} can be solved, one 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^+}$, and 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. A comparison of the kinetic and bond graph models is given in \autoref{fig:kinetic_BG_comp}. The bond graph model behaves similarly to the kinetic model. There is an initial difference between the models as the bond graph model moves to its steady state near $V=-120\si{mV}$, but the two models match very closely afterwards. At higher cycling velocities, the bond graph model has a slightly lower cycling velocity as the rapid equilibrium approximation starts to become less accurate. The CellML code describing the bond graph model is provided with this paper, and reproduces the curve in \autoref{fig:updated_original_comp}.

\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.} }
	\centering
	\begin{tabular}{c c}
		\toprule
		\textbf{Parameter} & \textbf{Value} \\ \midrule
			$k_1^+$ & $1423.2\si{s^{-1}}$\\ 
			$k_1^-$ & $225.9048\si{s^{-1}}$\\ 
			$k_2^+$ & $11564.8064\si{s^{-1}}$\\ 
			$k_2^-$ & $36355.3201\si{s^{-1}}$\\ 
			$k_3^+$ & $194.4506\si{s^{-1}}$\\ 
			$k_3^-$ & $281037.2758\si{mM^{-2}s^{-1}}$\\ 
			$k_4^+$ & $30629.8836\si{s^{-1}}$\\ 
			$k_4^-$ & $1.574\times 10^{6}\si{s^{-1}}$\\ 
			$K_{\text{d,Nai}}^0$ & $579.7295\si{mM}$\\ 
			$K_{\text{d,Nae}}^0$ & $0.034879\si{mM}$\\ 
			$K_{\text{d,Nai}}$ & $5.6399\si{mM}$\\ 
			$K_{\text{d,Nae}}$ & $10616.9377\si{mM}$\\ 
			$K_{\text{d,Ki}}$ & $16794.976\si{mM}$\\ 
			$K_{\text{d,Ke}}$ & $1.0817\si{mM}$\\ 
			$K_{\text{d,MgATP}}$ & $140.3709\si{mM}$\\ 
			$\Delta$ & -0.0550 \\
			Pump density & 1360.2624 \\ \bottomrule
	\end{tabular}
	\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. The reactions labelled 1,2,3 and 4 in the kinetic model correspond to the reactions R6, R7, R13 and R15 in the bond graph model respectively.}
	\centering
	\begin{tabular}{c c c}
		\toprule
		\textbf{Component} & \textbf{Parameter} & \textbf{Value} \\ \midrule
			R1 & $\kappa_1$ & $ 330.5462\si{fmol/s} $\\ 
			R2 & $\kappa_2$ & $ 132850.9145\si{fmol/s} $\\ 
			R3 & $\kappa_3$ & $ 200356.0223\si{fmol/s} $\\ 
			R4 & $\kappa_4$ & $ 2238785.3951\si{fmol/s} $\\ 
			R5 & $\kappa_5$ & $ 10787.9052\si{fmol/s} $\\ 
			R6 & $\kappa_6$ & $ 15.3533\si{fmol/s} $\\ 
			R7 & $\kappa_7$ & $ 2.3822\si{fmol/s} $\\ 
			R8 & $\kappa_8$ & $ 2.2855\si{fmol/s} $\\ 
			R9 & $\kappa_9$ & $ 1540.1349\si{fmol/s} $\\ 
			R10 & $\kappa_{10}$ & $ 259461.6507\si{fmol/s} $\\ 
			R11 & $\kappa_{11}$ & $ 172042.3334\si{fmol/s} $\\ 
			R12 & $\kappa_{12}$ & $ 6646440.3909\si{fmol/s} $\\ 
			R13 & $\kappa_{13}$ & $ 597.4136\si{fmol/s} $\\ 
			R14 & $\kappa_{14}$ & $ 70.9823\si{fmol/s} $\\ 
			R15 & $\kappa_{15}$ & $ 0.015489\si{fmol/s} $\\ 
			$\text{P}_1$ & $K_1$ & $101619537.2009\si{fmol^{-1}}$ \\ 
			$\text{P}_2$ & $K_2$ & $63209.8623\si{fmol^{-1}}$ \\ 
			$\text{P}_3$ & $K_3$ & $157.2724\si{fmol^{-1}}$ \\ 
			$\text{P}_4$ & $K_4$ & $14.0748\si{fmol^{-1}}$ \\ 
			$\text{P}_5$ & $K_5$ & $5.0384\si{fmol^{-1}}$ \\ 
			$\text{P}_6$ & $K_6$ & $92.6964\si{fmol^{-1}}$ \\ 
			$\text{P}_7$ & $K_7$ & $4854.5924\si{fmol^{-1}}$ \\ 
			$\text{P}_8$ & $K_8$ & $15260.9786\si{fmol^{-1}}$ \\ 
			$\text{P}_9$ & $K_9$ & $13787022.8009\si{fmol^{-1}}$ \\ 
			$\text{P}_{10}$ & $K_{10}$ & $20459.5509\si{fmol^{-1}}$ \\ 
			$\text{P}_{11}$ & $K_{11}$ & $121.4456\si{fmol^{-1}}$ \\ 
			$\text{P}_{12}$ & $K_{12}$ & $3.1436\si{fmol^{-1}}$ \\ 
			$\text{P}_{13}$ & $K_{13}$ & $0.32549\si{fmol^{-1}}$ \\ 
			$\text{P}_{14}$ & $K_{14}$ & $156.3283\si{fmol^{-1}}$ \\ 
			$\text{P}_{15}$ & $K_{15}$ & $1977546.8577\si{fmol^{-1}}$ \\ 
			$\text{K}_\text{i}^+$ & $K_\text{Ki}$ & $0.0012595\si{fmol^{-1}}$ \\ 
			$\text{K}_\text{e}^+$ & $K_\text{Ke}$ & $0.009236\si{fmol^{-1}}$ \\ 
			$\text{Na}_\text{i}^+$ & $K_\text{Nai}$ & $0.00083514\si{fmol^{-1}}$ \\ 
			$\text{Na}_\text{e}^+$ & $K_\text{Nae}$ & $0.0061242\si{fmol^{-1}}$ \\ 
			$\text{MgATP}$ & $K_\text{MgATP}$ & $2.3715\si{fmol^{-1}}$ \\ 
			$\text{MgADP}$ & $K_\text{MgADP}$ & $7.976 \times 10^{-5} \si{fmol^{-1}}$ \\ 
			$\text{P}_\text{i}$ & $K_\text{Pi}$ & $0.04565\si{fmol^{-1}}$ \\ 
			$\text{H}^+$ & $K_\text{H}$ & $0.04565\si{fmol^{-1}}$ \\ 
			Membrane capacitance & $C_m$ & $153400\si{fF}$ \\
			zF{\_}5 & $z_5$ & $-0.0550$ \\
			zF{\_}8 & $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=0.9\linewidth]{Figures/Terkildsen_NaK_abg}
	\caption{\textbf{Overall bond graph structure of the Terkildsen et al. model.}}
	\label{fig:Terkildsen_NaK}
\end{figure}

\begin{figure}
	\centering
	{\Large (a) \textbf{\textsf{P1{\_}6}}} \\
	\includegraphics[width=\linewidth]{Figures/P1_6_abg} \\[1cm]
	{\Large (b) \textbf{\textsf{P8{\_}13}}} \\
	\includegraphics[width=\linewidth]{Figures/P8_13_abg}\\[1cm]
	{\Large (c) \textbf{\textsf{P14{\_}15}}}\\
	\includegraphics[width=0.6\linewidth]{Figures/P14_15_abg}\\
	\caption{\textbf{Bond graph submodules in \autoref{fig:Terkildsen_NaK}.} (a) \textbf{\textsf{P1{\_}6}} submodule; (b) \textbf{\textsf{P8{\_}13}} submodule; (c) \textbf{\textsf{P14{\_}15}} submodule}
	\label{fig:BG_submodules}
\end{figure}

\end{document}