- Author:
- overleaf <overleaf@localhost>
- Date:
- 2017-11-02 00:40:20+00:00
- Desc:
- Responded to Kenneth's comments
- Permanent Source URI:
- https://models.physiomeproject.org/workspace/578/rawfile/bf2cb225936f91ae2151b278308613634a37ffa6/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}
\citet{terkildsen_balance_2007} developed a model of the Na$^+$/K$^+$ ATPase in cardiomyocytes to account for biophysics and thermodynamics of ion transport. 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 Na$^+$/K$^+$ ATPase model presented in \citet{terkildsen_balance_2007} model was based on that of \citet{smith_development_2004}, and fits the model to a wider range of data. It follows the methods proposed by \citet{smith_development_2004} to incorporate thermodynamic constraints and to simplify the model by using a lumping scheme. The Na$^+$/K$^+$ ATPase 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_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 the Na$^+$/K$^+$ ATPase model described in \citet{terkildsen_balance_2007} as the Terkildsen et al. model.
A key advantage of the Terkildsen et al. model over the earlier model by \citet{smith_development_2004} 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 was also activated by intracellular sodium, extracellular potassium, and MgATP. Data describing the activation of the Na$^+$/K$^+$ ATPase by 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 the same data set (\autoref{sec:modifications}). To verify the physical plausibility of the updated model, we developed 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 for 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 the cycle flux, with metabolite dependence incorporated in a manner that accounted for thermodynamic constraints.
We found three issues in the Terkildsen et al. 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}^0$ 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}^0$.
\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 expressions 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 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}
Because it was not possible to fix the above issues without significantly changing the kinetics of the model, we reparametersised 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} &
% * <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.
%
% ^.
\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 following 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 similar methods as Terkildsen et al. \citep{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 process 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 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 quality of fit is comparable to that in \citet{terkildsen_modelling_2006}. The fit to data for lower concentrations of extracellular sodium (\autoref{fig:fitting}(a)) was slightly poorer than the original model, although the unnormalised cycling velocities (\autoref{fig:fitting}(b)) appear to be more consistent with experimental data that suggest the saturated cycling velocity at positive membrane potentials is relatively insensitive to extracellular sodium \citep{nakao_[na]_1989}. The parameters of the updated model are given in \autoref{sec:parameters}.
% * <ktra014@aucklanduni.ac.nz> 2017-09-22T03:14:23.497Z:
%
% Add a sentence here comparing the fit of the updated model to the original model. Was the fit good just as good?
%
% ^ <panm@student.unimelb.edu.au> 2017-09-25T01:22:12.193Z.
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 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. We compared scaled versions of the updated and original 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 more 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).
% * <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 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 l}
\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 l}
\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}
