Computing Biomolecular System Steady-states
Peter Gawthrop \IEEEmembership{Senior Member
IEEE
Peter J. Gawthrop\footnote{Corresponding author. peter.gawthrop@unimelb.edu.au
A new approach to computing the equilibria and steady-states of biomolecular systems modelled by bond graphs is presented. The approach is illustrated using a model of a biomolecular cycle representing a membrane transporter and a model of the mitochondrial electron transport chain.
arxiv.tex
%%\RequirePackage{snapshot}
%%\documentclass[journal,a4paper]{IEEEtran}
%%\documentclass[draftcls,peerreview,11pt,a4paper]{IEEEtran}
\documentclass[a4paper.12pt]{article}
\title{Computing Biomolecular System Steady-states} 
%%\author{Peter Gawthrop \IEEEmembership{Senior Member, IEEE}}

\usepackage[noblocks]{authblk}
\author{Peter J. Gawthrop\footnote{Corresponding author. \textbf{peter.gawthrop@unimelb.edu.au}}}
\affil{
  Systems Biology Laboratory,
  Department of Biomedical Engineering,
  Melbourne School of Engineering,
  University of Melbourne,
  Victoria 3010, Australia.
   }

\usepackage{a4wide}

%Font
%%\usepackage{times}

%%Citations
\usepackage[numbers]{natbib}
\bibliographystyle{unsrtnat}

%%Maths
\usepackage{amsmath,amssymb,amscd,amstext,mathtools,extarrows,centernot}
\newcommand{\asinh}{\text{ asinh }}
%%\setcounter{MaxMatrixCols}{20}
\usepackage{chemformula}

%% Stoichiometric
\newcommand{\NN}{N}
\newcommand{\Nf}{{\NN^f}}
\newcommand{\Nr}{{\NN^r}}
\newcommand{\GG}{G}
%\newcommand{\KK}{K}
\newcommand{\Ncd}{\NN^{cd}}
\newcommand{\VV}{V}


\newcommand{\vv}{v_0}
\newcommand{\svv}{\vv^\star}
\newcommand{\svvc}{\vv^0}
\newcommand{\svveq}{\vv^{eq}}
\newcommand{\vp}{\vv^+}
\newcommand{\vm}{\vv^-}

\newcommand{\VN}{\phi_N}
\newcommand{\sPhi}{\Phi^\star}
\newcommand{\skappa}{\kappa^\star}
\newcommand{\skappac}{\kappa^0}
\newcommand{\skappaeq}{\kappa^{eq}}
%%\newcommand{\kappaeq}{\kappa^{eq}}
\newcommand{\kappabar}{\bar{\kappa}}

%% Electrical
\newcommand{\elec}{\text{e}^{-}}
\newcommand{\Vo}{\text{V}}

% \newcommand{\FV}{\text{V}^\prime}
% \newcommand{\FA}{\text{A}^\prime}
\newcommand{\FV}{\mathbb{V}}
\newcommand{\FA}{\mathbb{A}}
% \newcommand{\FV}{\text{\(\mathbb{V}\)}}
% \newcommand{\FA}{\text{\(\mathbb{A}\)}}

\newcommand{\tphi}{{\tilde{\phi}}}
\newcommand{\tPhi}{{\tilde{\Phi}}}

% Species
\usepackage{mhchem}
\newcommand{\pH}{\text{pH}}
\newcommand{\HATP}{\text{HATP}^{3-}}
\newcommand{\ADP}{\text{ADP}^{3-}}
\newcommand{\ATP}{\text{ATP}^{4-}}
\newcommand{\iP}{\text{Pi}}
\newcommand{\MgATP}{\text{MgATP}^{2-}}
\newcommand{\KATP}{\text{ATP}^{3-}}
\newcommand{\HPO}{\text{HPO}_4^{2-}}
\newcommand{\PI}{\text{P}_i^{2-}}
\newcommand{\POOO}{\text{PO}_3^{ -}}
\newcommand{\CT}{\text{C}_T}
\newcommand{\CoQ}{\text{CoQ}}
\newcommand{\Cox}{\text{C(ox)}^{3+}}
\newcommand{\Cred}{\text{C(red)}^{2+}}
\newcommand{\Feox}{\text{Fe}^{3+}}
\newcommand{\Fered}{\text{Fe}^{2+}}
\newcommand{\Fe}{\text{Fe}}
\newcommand{\Mnred}{\text{Mn}^{2+}}
\newcommand{\Mn}{\text{Mn}}
\newcommand{\Mnox}{\text{Mn}^{3+}}
\newcommand{\HHH}{\text{H}_2}
\newcommand{\HHO}{\text{H}_2\text{O}}
\newcommand{\HHHO}{\text{H}_3\text{O}^{+}}
\newcommand{\HHOO}{\text{H}_2\text{O}_2}
\newcommand{\HH}{\text{H}^{+}}
\newcommand{\Hi}{\text{H}_i^{+}}
\newcommand{\Hx}{\text{H}_x^{+}}
\newcommand{\KK}{\text{K}^{+}}
\newcommand{\Mg}{\text{Mg}^{2+}}
\newcommand{\NADH}{\text{NADH}}
\newcommand{\NAD}{\text{NAD}^{+}}
\newcommand{\OO}{\text{O}_2}
\newcommand{\sOO}{\text{O}^{\bullet -}_2}
\newcommand{\PP}{\text{P}}
\newcommand{\QH}{\text{QH}_2}
\newcommand{\QT}{\text{QH}_T}
\newcommand{\Q}{\text{Q}}
\newcommand{\ee}{\text{e}^{-}}

%Maths
\usepackage{amsmath,amssymb,amscd,amstext,mathtools,extarrows}
%%\numberwithin{equation}{section}

\newcommand{\hp}{{\cdot}}         %Hadamard or Schur product.
%%\newcommand{\hp}{\circ}         %Hadamard or Schur product.
\newcommand{\half}{\frac{1}{2}}  
\newcommand{\xt}{\tilde{x}}
\newcommand{\qt}{\tilde{q}}
\newcommand{\phit}{\tilde{\phi}}
\newcommand{\Phit}{\tilde{\Phi}}
\newcommand{\Phib}{{\bar{\Phi}}}
\newcommand{\Faraday}{\mathbf{F}}
\newcommand{\pmf}{\Delta p}
\newcommand{\emf}{\Delta \psi}

%%\newcommand{\Std}{\varnothing}
\newcommand{\Std}{\ominus}
\newcommand{\std}{\oslash}
\newcommand{\dx}{\dot{x}}
\newcommand{\dX}{\dot{X}}
\newcommand{\Gcm}{G_{cm}}
\newcommand{\phicm}{\phi_{cm}}
\newcommand{\lb}{\left (}
\newcommand{\rb}{\right )}
\newcommand{\reacul}[2]{
  {\; \xrightleftharpoons[#2]{#1} \;}
}


\newcommand{\reacu}[1]{
  \reacul{#1}{}
}

\newcommand{\reac}{
  \reacu{}
}

%Scientific units
\usepackage{siunitx}
\sisetup{per-mode = symbol}

%% Bond-Graph commands
%%\usepackage{BG}
 %% Bond-Graph commands
\newcommand{\BG}[1]{\text{\sffamily\textbf{#1}}}
\newcommand{\AmpE}{\BG{AE }}
\newcommand{\AF}{\BG{AF }}
\newcommand{\C}{\BG{C }}
\newcommand{\CS}{\BG{CS }}
\newcommand{\CSW}{\BG{CSW }}
\newcommand{\CSWw}{\BG{CSW}}
\newcommand{\De}{\BG{De }}
\newcommand{\Df}{\BG{Df }}
\newcommand{\I}{\BG{I }}
\newcommand{\INTF}{\BG{INTF }}
\newcommand{\ISW}{\BG{ISW}}
\newcommand{\R}{\BG{R }}
\newcommand{\Se}{\BG{Se }}
\newcommand{\See}{\BG{Se}}
\newcommand{\Sf}{\BG{Sf }}
\newcommand{\Sff}{\BG{Sf}}
\newcommand{\Sv}{\BG{Sv }}
\renewcommand{\SS}{\BG{SS }}
\newcommand{\one}{\BG{1 }}
\newcommand{\zero}{\BG{0 }}
\newcommand{\One}{\BG{1}}
\newcommand{\Zero}{\BG{0}}
\newcommand{\zeroo}{\BG{0 }}
\newcommand{\TF}{\BG{TF }}
\newcommand{\FF}{\BG{F }}
\newcommand{\GY}{\BG{GY }}

\renewcommand{\Re}{\BG{Re }}
\newcommand{\mRe}{\BG{mRe }}
\newcommand{\sRe}{\BG{sRe }}
\newcommand{\iRe}{\BG{iRe }}
\newcommand{\Rr}{\BG{Rr }}
\newcommand{\Rrr}{\BG{Rr2 }}

%% Bond-Graph components with label
\newcommand{\BGL}[2]{$\BG{#1}$:$\mathbf{#2}$} %Generic

\newcommand{\BC}[1]{\BGL{C}{#1}}
\newcommand{\BCS}[1]{\BGL{CS}{#1}}
\newcommand{\BI}[1]{\BGL{I}{#1}}
\newcommand{\BR}[1]{\BGL{R}{#1}}

\newcommand{\BSS}[1]{\BGL{SS}{#1}}
\newcommand{\BSe}[1]{\BGL{Se}{#1}}
\newcommand{\BSf}[1]{\BGL{Sf}{#1}}
\newcommand{\BSv}[1]{\BGL{Sv}{#1}}
\newcommand{\BDe}[1]{\BGL{De}{#1}}
\newcommand{\BDf}[1]{\BGL{Df}{#1}}

\newcommand{\BTF}[1]{\BGL{TF}{#1}}
\newcommand{\BGLY}[1]{\BGL{GY}{#1}}

\newcommand{\BRe}[1]{\BGL{Re}{#1}}
\newcommand{\BmRe}[1]{\BGL{mRe}{#1}}
\newcommand{\BsRe}[1]{\BGL{sRe}{#1}}
\newcommand{\BiRe}[1]{\BGL{iRe}{#1}}



%Code
\usepackage{listings}

%% Figs
\usepackage{graphicx,subfigure,fancybox,color}
\newcommand{\FigWidth}{0.75}
\newcommand{\Fig}[2]{
  \includegraphics[width=#2\linewidth]{#1.pdf}
   \label{subfig:#1}
}

\newcommand{\SubFig}[3]{
  \subfigure[#2]{
    \includegraphics[width=#3\linewidth]{#1.pdf}
    \label{subfig:#1}
  }
}

\newcommand{\SubFigH}[3]{
  \subfigure[#2]{
    \includegraphics[height=#3\linewidth]{#1.pdf}
    \label{subfig:#1}
  }
}

%%\usepackage[pagebackref]{hyperref}
\usepackage{hyperref}
\hypersetup{colorlinks=true,citecolor=blue,linkcolor=blue,urlcolor=green}
\usepackage{url}
\usepackage{doi}


\begin{document}
\maketitle
\begin{abstract}
  A new approach to computing the equilibria and steady-states of
  biomolecular systems modelled by bond graphs is presented. The
  approach is illustrated using a  model of a biomolecular cycle
  representing a membrane transporter and a model of the mitochondrial
  electron transport chain.
\end{abstract}

% \newpage
% \tableofcontents
% \newpage

% \section*{TBD}
% \begin{enumerate}
% \item Use $\frac{\phi}{V_N}$ as in IET paper.
% \item Sort out molar fac v mols.
% \item Test stoich\_equilibrium on Ox examples  (already tesed on ABE examples).
% \item New stoich\_steady to replace old stoich\_ss.
% \end{enumerate}
\section{Introduction}
\label{sec:introduction}
%% Modelling
There are many ways of deriving the dynamical equations describing the
time course of flows and concentrations in biomolecular systems
\citep{Bea11,Bea12}.  As physical systems, biomolecular systems are
subject to the laws of thermodynamics \citep{BeaQia10,AtkPau11} and
number of such \emph{energy-based} approaches have been developed
including the bond graph of \citet{OstPerKat71,OstPerKat73} further
developed by \citet{GreCel12} and
\citet{GawCra14,GawCra16,GawCra17}. The bond graph approach is
particularly appropriate for systems which combine different physical
domains and thus are particularly appropriate for electrogenic
biomolecular systems \citep{GawSieKam17,Gaw17a}. In common with other
approaches, the bond graph approach can be used to derive ordinary
differential equations (ode) describing the time evolution of the
chemical species, electrical charges and chemical and electrical
flows.

%% Bond Graph  Material
The basic ideas of bond graphs can be found in textbooks
\citep{Bor10,KarMarRos12} and a tutorial for engineers \citep{GawBev07}.
%
Appendices A \& B of a recent paper \citep{GawCra17} provide a
tutorial introduction to bond graphs in the context of biomolecular
systems; in particular, the stoichiometric matrix describing the
biomolecular system can be automatically derived from the bond graph. In
this paper, both the bond graph and relevant stoichiometric matrices
are stated but the derivation is omitted to save space and to focus
attention on the novel ideas of the paper.

%% Equilibria
In some cases, it is of interest to find values corresponding to a
situation where the amount of each chemical species and electrical
charge is constant; this is termed a \emph{steady-state} of the
system. System steady-states correspond to constant chemical and
electrical flows; in the special case the such flows are zero,
the system is said to be in \emph{equilibrium}.
%
For example, study of the energetic behaviour of metabolic
biomolecular systems at constant flow leads to conclusions about
efficient enzyme utilisation \citep{ParRubXu16}.

%
A number of approaches to finding the steady-state of biomolecular
networks have been suggested including simulation, flux balance
analysis~\citep{OrtThiPal10} and the diagram methods of
\citet{KinAlt56} and \citet{Hil89}.
 
%%The simulation approach.
The simulation approach is exemplified by \citet{BazBeaVin16} who
state that ``For the steady state analysis, the model was simulated
until the absolute values of the state variable derivatives were less
than $10^{-10}$ ...''.
%
A similar approach is to simulate starting from equilibrium with
slowly-varying chemostats followed by an algebraic solver for each
simulation time point starting with the simulated state
\citet{GawCra17}. Such approaches can be numerically challenging and time
consuming.
%
%% Flux Balance Analysis
Flux balance analysis (FBA) \citep{OrtThiPal10} is based on the system
stoichiometric matrix but  ignores reaction kinetics.
Essentially, the system
stoichiometric matrix constrains all steady-state flows
$\VV$ to lie within the right null subspace of the chemodynamic
stoichiometric matrix $\Ncd$. 
%
These constraints are not enough to deduce the steady-state flows
$\VV$. Specific values for $\VV$ are found by adding extra constraints
together with a quadratic (in the flows) cost function and applying
linear programming.
%
Energy considerations are not directly addressed, but can be added on
\citep{HenBroHat07}.
%
The diagram method of \citet{KinAlt56} and \citet{Hil89} can be used
to generate explicit rate equations for reaction networks using, for
example, the software described by \citet{QiDasHan09}.

A general bond graph approach to steady state analysis using causality
concepts was introduced by \citet{Bre84b} and extended by
\citet{BidMarSca07} and \citet{GonGal11}. In contrast, the approach of
this paper is specifically tailored to the special properties of
biomolecular system bond graphs: no inertial components, conserved
moieties leading to initial condition dependent steady-states and
particular non-linearities involving logarithms and exponentials.
%
The special bond graph junction structure of biochemical systems is
noted elsewhere \citep{Gaw17}. In particular, in a similar fashion to
\citep{Bre84b}, the application of derivative causality reveals key
stoichiometric properties such as conserved moieties and flux
pathways. 

The approach to steady-state analysis proposed here is based on two
ideas: using the \emph{Faraday equivalent potential} \citep{Gaw17a} to
unify chemical and electrical potential and using the positive-pathway
approach \citep{GawCra17} to represent the steady-state flow pathways.

%% New method
% The approach suggested herein is summarised as follows:
% \begin{enumerate}
% \item Use the positive-pathway approach of \citet{GawCra17} to
%   identify pathways: typically there will be few pathways even for
%   complex biomolecular networks.
% \item Choose the pathway flows $v_p$.
% \item Use an iterative approach to deduce reaction affinities $\Phi$ from
%   flows.
% \item Compute the potentials $\phi$ corresponding to $\Phi$ constrained
%   by chemostats and conserved moieties.
% \end{enumerate}
The approach suggested here is like the FBA approach insofar as it
uses pathway analysis of steady-state flows, and like the simulation
approach in that it uses reaction kinetics. When based on bond graph
models, the approach ensures thermodynamic compliance.

\S~\ref{sec:background} gives the equations required in the rest of
the paper.
%
\S~\ref{sec:reaction-kinetics} shows how conventional mass-action
kinetics, expressed in terms of forward $\Phi^f$ and reverse $\Phi^r$
reaction affinities, can be re-expressed in terms of reaction affinity
$\Phi = \Phi^f - \Phi^r$ and summed reaction affinity
$\sPhi = \Phi^f + \Phi^r$. This involves the hyperbolic functions
$\sinh$ and $\tanh$.
%
\S~\ref{sec:equilibria} shows how the equilibrium values of the  free
potentials can be deduced from the fixed potentials.
%
\S~\ref{sec:steady-states} contains the main results of the paper:
how the steady-state values of the free potentials can be deduced from
the fixed potentials for a given pathway flow.
%
\S~\ref{sec:exampl-biom-cycle} uses the biomolecular cycle of
\citet{Hil89} to illustrate the approach.
%
\S~\ref{sec:exampl-oxid-phosph} uses the Mitochondrial Electron Transport
Chain as a more detailed example.
%
\S~\ref{sec:conclusion} concludes the paper and suggest future research directions.

\section{Background}
\label{sec:background}
This section contains the background material required in the rest of
the paper. More information is to found in a companion
paper~\citep{Gaw17a}.

The Faraday-equivalent potential~\citep{Gaw17a} of a substance A is
given in terms of the potential $\phi_A^\std$ (with units of
Volts~\si{V}) when the amount of the substance is $x_A^\std\si{mol}$
as:
\begin{xalignat}{2}
  \phi_A &= \phi_A^\std + \phi_N \ln \frac{x_A}{x^\std_A} \label{eq:phi_A}&
\text{where } 
  \phi_N &= \frac{RT}{F} %\label{eq:V_N}
\end{xalignat}
$R = \SI{8.314}{JK^{-1}mol^{-1}}$ is the universal gas constant,
$T~\si{K}$ is the absolute temperature and
$F \approx \SI{96485}{C.mol^{-1}}$ is Faraday's constant.

%%\paragraph{Remarks}
\begin{enumerate}
% \item The Faraday-equivalent  potential $\phi$ has units of volts
%   (\si{V}) and
%   the quantity $x_A$ has units of moles (\si{mol}).
\item For a given temperature, $\phi_N$ is a constant given by
  \eqref{eq:phi_A} and, at $T = \SI{300}{C}$ has a value of about
  $\SI{26}{mV}$.
\item Equation \eqref{eq:phi_A} has two parameters:
  \begin{enumerate}
  \item $\phi_A^\std$ the potential of the substance in the relevant
    compartment at nominal conditions and
  \item $x^\std_A$, the amount of the substance in the relevant
    compartment at the nominal
    conditions defining the nominal potential $\phi_A^\std$.
  \end{enumerate}
  % Thus the substance A appearing in two different compartments would
  % have the same nominal potential $\phi_A^\std$ but different
  % nominal amounts $x^\std_A$.
\item Given the two parameters $\phi_A^\std$ and $x^\std_A$, and the
  constant $\phi_N$, the amount $x_A$ may be deduced from $\phi_A$ by
  inverting Equation \eqref{eq:phi_A}:
  \begin{equation}
    x_A = x^\std_A \exp \frac{\phi_A - \phi_A^\std}{\phi_N} \label{eq:x_a}
  \end{equation}
\item The potential $\phi_A^\std$ of
  substance A at \emph{nominal} conditions can be derived from the
  potential $\phi_A^\Std$ at \emph{standard} conditions using the
  formula~\citep{AtkPau11}:
  \begin{equation}
    \phi_A^\std = \phi_A^\Std + \phi_N \ln \frac{x^\std_A}{x^\Std_A} 
  \end{equation}
Standard chermical potentials are listed for common substances by
\citet{AtkPau11}.
\item Equation \eqref{eq:phi_A} can be rewritten as:
  \begin{xalignat}{2}
    \phi_A &=  \phi_N \ln K_A x_A \label{eq:phi_A_K}&
    \text{where }
    K_A &= \frac{\exp \frac{\phi_A^\std}{\phi_N}}{x^\std_A}%\label{eq:K}
  \end{xalignat}
\end{enumerate}

The forward $\Phi^f$ and reverse $\Phi^r$ reaction affinities can be
deduced from the Faraday-equivalent  potentials $\phi$ of the
substances using the forward $\Nf$ and reverse $\Nr$ stoichiometric
matrices as:
\begin{xalignat}{2}
  \Phi^f &= \Nf^T \phi \label{eq:Phi^f}&
  \Phi^r &= \Nr^T \phi %\label{eq:Phi^r}
\end{xalignat}
Furthermore, the net reaction affinity $\Phi$ is defined as
\begin{align}
  \Phi &= \Phi^f - \Phi^r \label{eq:Phi_def}\\
  \text{ and so }
  \Phi &= -\NN^T \phi \label{eq:Phi}
 \text{ where } 
  \NN = \Nr - \Nf
\end{align}
and $N$ is the system \emph{stoichiometric matrix}.
%
$N$, $\Nf$ and $\Nr$ can be automatically deduced from the system bond
graph \citep{GawCra14}.

Mass-action reaction kinetics can be written as:
\begin{align}
  v &= \kappa \lb \vp - \vm \rb\label{eq:MA_0}\\
  \text{where }
  \vp &= \exp \frac{\Phi^f}{\VN} \label{eq:vp}\\
  \text{and }
  \vm &= \exp \frac{\Phi^r}{\VN} \label{eq:vm}
\end{align}
where $\kappa$ is a constant and $\Phi^f$ and $\Phi^r$ are the forward
and reverse reaction affinities. 

From equations \eqref{eq:MA_0}~--~\eqref{eq:vm}, zero reaction flow $v=0$ implies that
\begin{xalignat}{2}
  \vp &= \vm & %\label{eq:vp_0}\\
\text{and }
  \Phi &= \Phi^f - \Phi^r =0 \label{eq:Phi_0}
\end{xalignat}


As discussed in the literature, the system ode can be written as:
\begin{equation}
  \dX = N V
\end{equation}
where $\dX$ is the rate of change of the vector of chemical species
and electrical charge, $N$ is the system stoichiometric matrix
and
$V$ is the vector of chemical and
electrical flows.  

The notion of a \emph{chemostat}~\citep{PolEsp14,GawCra16}, the
chemical analogue of a voltage source, is useful when considering the
steady-states of open biomolecular systems.
%
In particular, as discussed by \citet{GawCra16}, in the presence of
chemostats:
\begin{equation}\label{eq:dX_cs}
  \dX = \Ncd V
\end{equation}
where $\Ncd$ is $N$ with the rows corresponding to the chemostats set
to zero.  In either case, the reaction potentials $\Phi$ are given by
Equation \eqref{eq:Phi}.

The left null space matrix $G$ of $\Ncd$ has the property that:
\begin{align}
  G\Ncd &= 0\\
  \text{hence }
  G\dX &= G\Ncd V = 0
\end{align}
The $n_{cm}$ rows of $G$ thus define the $n_{cm}$ \emph{conserved
  moieties} given by:
\begin{align}\label{eq:cm}
  GX = GX_0
\end{align}
where $X_0$ is the initial state.

In this paper, the \emph{order}  of the $i$th conserved moiety is
defined as the the number of non zero entries in the $i$th row of $G$
minus $1$. Thus zero-order moieties correspond to chemostats and first
order moieties to two substances whose net amount is constant.

% In the zero-order case, all entries of $G$ can be set to unity by
% suitably redefining $X^{cm}$. This will be assumed in the sequel.

\section{Reaction Kinetics}
\label{sec:reaction-kinetics}
% Mass-action reaction kinetics of eqn. \eqref{eq:MA_0} are such that
% the flow is zero if $\Phi^f=\Phi^r$ or $\Phi^f-\Phi^r=0$.
%
The mass-action kinetics of eqns.~\eqref{eq:MA_0}~--~\eqref{eq:vm} are expressed in
terms of the forward and reverse reaction affinities $\Phi^f$ and
$\Phi^r$. One property of eqns.~\eqref{eq:MA_0}~--~\eqref{eq:vm} is that flow $v$ is
zero if $\Phi^f=\Phi^r$ or $\Phi=\Phi^f-\Phi^r=0$ where $\Phi$ is the
(net) reaction affinity.

To obtain further insight into mass-action reaction kinetics it is
therefore convenient to rewrite eqn.~\eqref{eq:MA_0} explicitly in terms of the
reaction affinity $\Phi$ and \emph{summed reaction affinity} $\sPhi$ where:
\begin{xalignat}{2}
  \Phi &= \Phi^f - \Phi^r&
  \text{and }
  \sPhi &= \Phi^f + \Phi^r\label{eq:sPhi}
\end{xalignat}
These definitions imply that 
\begin{xalignat}{2}
  \Phi^f &= \frac{\sPhi + \Phi}{2}&
  \text{and }
  \Phi^r &= \frac{\sPhi - \Phi}{2}
\end{xalignat}
and hence that
\begin{xalignat}{3}
  \vp &= \skappa \exp \frac{\Phi}{2\VN} &
  \text{and }
  \vm &=  \skappa \exp \frac{-\Phi}{2\VN}\\
  \text{where }
  \skappa &= \exp \frac{\sPhi}{2\VN}\label{eq:skappa}
\end{xalignat}
Introducing the difference $\vv$ and sum $\svv$ of $\vp$ \eqref{eq:vp}
and $\vm$ \eqref{eq:vm}:
\begin{align}
  \vv &= \vp - \vm = 2\skappa \sinh \frac{\Phi}{2\VN}\label{eq:v_0}\\
  \svv &= \vp + \vm = 2\skappa \cosh \frac{\Phi}{2\VN}\label{eq:sv_0}
\end{align}
where $\sinh$ and $\cosh$ are the hyperbolic $\sin$ and $\cos$
functions respectively.

Using these equations, the mass-action kinetics of eqn.~\eqref{eq:MA_0} can be
rewritten in two alternative forms:
\begin{xalignat}{1}
  v &= 2\kappa \skappa \sinh \frac{\Phi}{2\VN} \label{eq:MA_sinh}\\
%%&  \text{where } \skappa &= 2\exp \frac{\Phi^f + \Phi^r}{2\VN}\\
  v &= \kappa \svv \tanh \frac{\Phi}{2\VN}   \label{eq:MA_tanh}
\\
%%&
  \text{where } \skappa &= \exp \frac{\Phi^f + \Phi^r}{2\VN}\label{eq:skappa_1}\\
  \text{and } \svv &= \exp \frac{\Phi^f}{\VN} + \exp \frac{\Phi^r}{\VN}\label{eq:svv}
\end{xalignat}
$\tanh$ is the hyperbolic $\tan$ function, 
%
\eqref{eq:skappa_1} follows from \eqref{eq:sPhi}\&\eqref{eq:skappa}
%
and \eqref{eq:svv} follows from Equations
\eqref{eq:vp}\&\eqref{eq:vm}.
%
Equation \eqref{eq:MA_sinh} is explicitly in terms of the reaction
affinity $\Phi$ \eqref{eq:Phi} and, via $\skappa$ \eqref{eq:skappa},
$\sPhi$ \eqref{eq:sPhi} and via \eqref{eq:Phi^f}, the potentials
$\phi$.
%
Equation \eqref{eq:MA_tanh} is also explicitly in terms of  the reaction
affinity $\Phi$ \eqref{eq:Phi} but, unlike Equation
\eqref{eq:MA_sinh}, depends on $\svv$ rather than $\skappa$.

Define $\skappaeq$ and $\svveq$ as the equilibrium values of
$\skappa$ and $\svv$ respectively.
%
At equilibrium, $\Phi^f=\Phi^r$ and thus, from Equations
\eqref{eq:skappa_1} and \eqref{eq:svv} it follows that
\begin{equation}
  \svveq = 2\skappaeq 
\end{equation}

Equations \eqref{eq:MA_sinh} and
\eqref{eq:MA_tanh} can be rewritten as:
\begin{xalignat}{2}
  v &=  2\kappabar \rho_s\sinh \frac{\Phi}{2\VN} &
  \text{where } 
   \rho_s &= \frac{\skappa}{\skappaeq}\label{eq:MA_sinh_0}\\
  v &= 2\kappabar \rho_t\tanh \frac{\Phi}{2\VN} &
  \text{where } 
  \rho_t &= \frac{\svv}{\svveq}
  \label{eq:MA_tanh_0}
\end{xalignat}
and $\kappabar = \kappa \skappaeq = \kappa \frac{\svveq}{2}$.
%
Near equilibrium, $\Phi \approx 0$ and these two equations can each be
approximated by functions of $\Phi$ only as:
\begin{align}
  v &\approx  2\kappabar \sinh \frac{\Phi}{2\VN}\label{eq:MA_sinh_a}\\
  v &\approx  2\kappabar \tanh \frac{\Phi}{2\VN}\label{eq:MA_tanh_a}
\end{align}
They may be further approximated by the linear equation:
\begin{align}
  v &\approx  \kappabar \frac{\Phi}{\VN} 
\end{align}

\subsection{Example: A Basic Electrogenic Network}
\label{sec:exampl-basic-electr}
\begin{figure}[htbp]
  \centering
  \Fig{ABE_cbg}{0.6}
  %%\SubFig{ABE_phi_E}{$\phi_E$}{0.45}
  % \SubFig{ABE_kappa}{$\kappa$}{0.45}
  % \SubFig{ABE_crite}{$\Delta$}{0.45}
  \caption[A Basic Electrogenic Network]{A Basic Electrogenic
    Network. The system bond graph has two \C components: \BC{A} and
    \BC{B} representing the chemical species A and B respectively and
    a third \C component \BC{E} representing electrical trans-membrane
  charge. The \Re component represents the reaction and the \zero and
  \one junction connect components appropriately. The bonds have been
  annotated to show the corresponding Faraday-equivalent potential
  (effort variable) and Faraday-equivalent current (flow
  variable). Thus the bond impinging on the \BRe{r} component has the
  effort variable $\Phi^f=\phi_A-\phi_E$ and flow variable $v$.
  % 
  The constitutive relations for the \C components are given by
  \eqref{eq:phi_A} and for the \Re component by
  \eqref{eq:MA_0}~--~\eqref{eq:vm}.
}
\label{fig:ABE_BG}
\end{figure}
Figure \ref{fig:ABE_BG} gives the bond graph  of a basic
electrogenic network where the trans-membrane reaction:
\begin{equation}
  A \reacu{r} B
\end{equation}
accumulates charge on the capacitor \BC{E}. 
For illustration, it is assumed that:
\begin{equation}
  \phi_A^\std = \phi_B^\std = 0
\end{equation}

The state vector is chosen as
\begin{align}
  X &=
      \begin{pmatrix}
        x_A&x_B&x_E
      \end{pmatrix}^T
\end{align}
where $x_A$, $x_B$ and $x_E$ are the amounts of $A$, $B$ and charge
respectively.  The stoichiometric matrices $N$, $\Nf$ and $\Nr$ are
automatically derived from the bond graph of Fig. \ref{fig:ABE_BG} and
are given by:
\begin{xalignat}{3}
  N &=
  \begin{pmatrix}
    -1\\1\\2
  \end{pmatrix}&
  \Nf &=
  \begin{pmatrix}
    1\\0\\-1
  \end{pmatrix}&
  \Nr &=
  \begin{pmatrix}
    0\\1\\1
  \end{pmatrix}
\end{xalignat}
Hence the corresponding reaction affinities are given by Equations
\eqref{eq:Phi^f}~--~\eqref{eq:Phi}
as
\begin{xalignat}{2}
  \Phi^f &= \phi_A - \phi_E&
  \Phi^r &= \phi_B + \phi_E\label{eq:Phir_example}
\end{xalignat}
and the corresponding difference and summed reaction affinities are:
\begin{align}
  \Phi &= \phi_A - \phi_B - 2\phi_E&
  \sPhi &= \phi_A + \phi_B \label{eq:Phi_example}
\end{align}


\begin{figure}[htbp]
  \centering
  %%\SubFig{xx_CM}{X: case 1}{0.45}
  \SubFig{CR_CM}{$x_A+x_B=10.2$}{\FigWidth}
  %%\SubFig{xx_noCM}{X: case 2}{0.45}
  \SubFig{CR_noCM}{$x_B=5$}{\FigWidth}
  \caption[Non-electrogenic case.]{Non-electrogenic case. 
    Flow $v$ computed exactly using Octave \citep{EatBatHau15} and using the $\sinh$ and  $\tanh$
    approximations.
    % 
    (a)  $x_A+x_B=10.2$. The $\tanh$ approximation is exact in this case.
    % 
%    (c) Variation of states $x_A$ and $x_B$ when $x_B=5$.
    % 
    (b)  $x_B=5$. Neither approximation is exact in this case.
  }
\label{fig:RK-non-electrogenic}
\end{figure}
Consider first the non-electrogenic case where $\phi_E=0$. Then using
the formula \eqref{eq:phi_A} and eqn. \eqref{eq:sv_0} for chemical
potential:
\begin{align}
  \svv &= \exp \frac{\Phi^f}{\VN} + \exp \frac{\Phi^r}{\VN}
       = \exp \frac{\phi_A}{\VN} + \exp \frac{\phi_B}{\VN}\notag\\
       &= K_A x_A + K_B x_B
\end{align}
For the purposes of illustration, further assume that $K_A=K_B=1$ thus
\begin{equation}
  \svv = x_A + x_B
\end{equation}
As $A$ and $B$ form a conserved moiety,  $\svv$ is constant and thus
the $\tanh$ approximation \eqref{eq:MA_tanh_a} is exact.  Figure
\ref{subfig:CR_CM} shows the exact flow $v$ together with the
approximations \eqref{eq:MA_sinh_a} and \eqref{eq:MA_tanh_a} for this
case.

In contrast, Figure \ref{subfig:CR_noCM} shows the case where $x_b=0.5$
and $x_a$ varies. There is \emph{no} conserved moiety:
$x_A+x_B \ne constant$.  In this case neither approximation is exact,
though the $\tanh$ approximation is better than the $\sinh$
approximation.

\begin{figure}[htbp]
  \centering
  \Fig{ABE_phi_E}{\FigWidth}
  %%\SubFig{ABE_kappa}{$\dfrac{\skappa}{\skappaeq}$}{0.45}
  \caption{Electrogenic case with fixed concentrations. 
    % 
    Using Octave \citep{EatBatHau15}, the capacitor potential $\phi_E$
    plotted against reaction flow when $\phi_A=\SI{500}{mV}$ and
    $\phi_B=\SI{100}{mV}$. 
    % 
    The reaction constant \eqref{eq:MA_sinh_0} is chosen as 
    % 
    $\kappa_s = 0.01$.
    % 
    At $v=0$, $\phi_E = \frac{\phi_A - \phi_B}{2}=\SI{200}{mV}$; as
    $v$ increases, $\phi_E$ decreases as the inverse hyperbolic sin of
    \eqref{eq:phi_E_form}.}
\label{fig:RK-electrogenic}
\end{figure}
Now consider the electrogenic, but constant concentration, case where
$\phi_E$ varies but $\phi_A$ and $\phi_B$ are constant. It follows
from Equation \eqref{eq:Phi_example} that $\skappa$ is constant and
thus the $\sinh$ approximation \eqref{eq:MA_sinh_a} is exact.
Combining equations \eqref{eq:MA_sinh_a} and \eqref{eq:Phi_example},
it follows that:
\begin{align}\label{eq:phi_E_form}
  \phi_E &= \frac{\phi_A - \phi_B}{2} - \phi_N \asinh \frac{v}{2 \kappabar}
\end{align}
This is plotted in Figure \ref{fig:RK-electrogenic} for a particular
choice of parameters. This result is qualitatively similar to the
experimentally obtained results of \citet[Figure 2B]{BazBeaVin16} for
the Mitochondrial Electron Transport Chain. This is explored further
in \S~\ref{sec:exampl-oxid-phosph}.

\section{Equilibria}
\label{sec:equilibria}
In typical biomolecular systems, the reaction kinetics of
\S~\ref{sec:reaction-kinetics} are  embedded in large networks. This
section considers the equilibria of such networks.
%
In this paper, a system is said to be in \emph{equilibrium} if efforts
(eg chemical potentials and voltages) are such that all flows are
zero,
\begin{equation}\label{eq:equilibrium}
  \VV = 0
\end{equation}
which implies from the arguments leading to equation \eqref{eq:Phi_0}
that the reaction affinities $\Phi=0$.

As discussed in \S~\ref{sec:background}, the affinities $\Phi$ are
given by  equation \eqref{eq:Phi} ($\Phi = -N^T \phi$).
% \begin{equation}\label{eq:Phiphi}
%   \Phi = -N^T \phi
% \end{equation}
Equation \eqref{eq:Phi} is linear and therefore finding the potentials
$\phi$ to give zero affinities $\Phi=0$ is a linear problem.  However,
the $n_V \times n_X$ matrix $N^T$ is usually not full rank as usually
$n_V < n_X$ and therefore \eqref{eq:Phi} cannot be used directly to
deduce the potentials $\phi$ from reaction affinities $\Phi$.

However, the conserved moieties of Equation \eqref{eq:cm} give a
further set of constraints on the potentials $\phi$. In particular,
using Equation \eqref{eq:x_a}, Equation \eqref{eq:cm} becomes:
\begin{align}
   G \lb X^\std \hp  \exp  \frac{\phi -  \phi^\std}{\phi_N} \rb &= GX_0
\label{eq:GX_phi}
\end{align}
where $\hp$ denotes elementwise multiplication.
Equation \eqref{eq:GX_phi} is \emph{nonlinear} in $\phi$. 


\subsection{Conserved Moieties with order zero}
In the special case that all conserved moieties are order $0$,
Equation \eqref{eq:GX_phi} can be rewritten as the linear equation.
\begin{align}
  G \phi &= \phi^{cm}\label{eq:Gphi}\\
  \text{where }
  \phi^{cm} & = G \phi_0
  \text{ and }
  \phi_0 = \phi^\std + \phi_N \ln \frac{X_0}{X^\std}
\end{align}
In this case Equations \eqref{eq:Phi} and  \eqref{eq:Gphi} can be
combined as:
\begin{equation}\label{eq:phi_cm_0}
  \begin{pmatrix}
    -N^T\\G
  \end{pmatrix}
  \phi =
  \begin{pmatrix}
    \Phi\\\phi^{cm}
  \end{pmatrix}
\end{equation}
Once again, equilibria are defined by setting $\Phi=0$ and the
solution of Equation \eqref{eq:phi_cm_0} is a linear problem. 

% The vector of potentials $\phi$ is partitioned into those associated
% with chemostats, and therefore fixed, and those that are free to vary:
% \begin{equation}
%   \phi =
%   \begin{pmatrix}
%     \phi_v \\
%     \hline
%     \phi_{cs} 
%   \end{pmatrix}
% \end{equation}
% Similarly, the stoichiometric matrix $\NN$ is partitioned conformabley
% as:
% \begin{equation}
%   \NN =
%   \begin{pmatrix}
%     \NN_v\\     
%     \hline 
%     \NN_{cs}
%   \end{pmatrix}
% \end{equation}
% Hence Equation \eqref{eq:Phi} can be rewritten as:
% \begin{align}
%   \Phi &= -\lb \NN_v^T \phi_v +  \NN_{cs}^T \phi_{cs}   \rb\\
% \text{or }
%   \NN_v^T \phi_v &= -\lb \NN_{cs}^T \phi_{cs} + \Phi \rb \label{eq:phi_cm_0}
% \end{align}
% In the case of equilibrium, $\Phi=0$ and
% \begin{align}
%   \NN_v^T \phi_v &= - \NN_{cs}^T \phi_{cs}  \label{eq:phi_v_eq}
% \end{align}
% Equation \eqref{eq:phi_v_eq} has not, a solution according to the usual rules of
% linear algebra. For example, if $\NN_v$ is square and of full rank:
% \begin{equation}\label{eq:phi_v_sq}
%  \phi_v = - \lb\NN_v^{T}\rb^{-1}\NN_{cs}^T \phi_{cs}
% \end{equation}

\subsection{Conserved Moieties with order greater than zero}
\label{sec:conserved-moieties}
The order-zero case gives the linear set of equations
\eqref{eq:phi_cm_0}. This section shows how the non-linear equations
arising when the conserved moiety is not zero-order may be solved.

Define $G_{cm}$ as the matrix containing the non-zero-order rows of $G$.
\begin{enumerate}
\item The non-zero-order moieties are converted to order-zero by deleting
  all but one element of the rows of $G_{cm}$.
\item A $\phicm$ for the new chemostats is chosen and the linear problem
  \eqref{eq:phi_cm_0} solved. 
\item Compute $\Gcm X$ and evaluate discrepancy with respect to $X^{cm}$.
\item Choose new $\phicm$ to reduce discrepancy and repeat from step 2.
\end{enumerate}
This is conveniently automated using \texttt{fsolve} in
Octave~\citep{EatBatHau15}.

\subsection{Example: A Basic Electrogenic Network
  (continued)}\label{sec:exampl-basic-electr-1}

\begin{figure}[htbp]
  \centering
  \SubFig{eq_phi_nocm}{$\phi_A$, $\phi_B$}{\FigWidth}
  \SubFig{eq_x_nocm}{$X$}{\FigWidth}
  \caption{Electrogenic system: Equilibria, order-zero moieties. A and
  E are chemostats; B is free to vary.
  % 
  (a) The potentials $\phi_A$ and $\phi_B$ when B is a chemostat with
  $\phi_B=0$ and the electrical potential $\phi_E$ varies.
  % 
  (b) The corresponding states $x_A$ and $x_B$.
}
\label{fig:nocm}
\end{figure}
Consider first the case where both A and E are chemostats. Then $\Ncd$
\eqref{eq:dX_cs} is given by $N$ with the first and third rows
deleted:
\begin{align}
 \Ncd &=
  \begin{pmatrix}
    0&1&0
  \end{pmatrix}^T
\end{align}
One choice of $G$ is
\begin{align}
  G &=
  \begin{pmatrix}
    1&0&0\\
    0&0&1
  \end{pmatrix}^T
\end{align}
As there is only one non-zero entry in each row of $G$,  the linear
equation \eqref{eq:phi_cm_0} can be used with:
\begin{xalignat}{2}
  \begin{pmatrix}
    -N^T\\G
  \end{pmatrix} &=
  \begin{pmatrix}
    -1&1&2\\
    1&0&0\\
    0&0&1
  \end{pmatrix};&
  \begin{pmatrix}
    \Phi\\\phi^{cm} 
  \end{pmatrix} &=
                  \begin{pmatrix}
                  0\\ \phi_A \\ \phi_E  
                  \end{pmatrix}
\end{xalignat}
Figure \ref{fig:nocm} shows the solution of equation
\eqref{eq:phi_cm_0} for $\phi_A$, $\phi_B$ plotted against $\phi_E$
when the potential chemostat A is fixed at $\phi_A=0$ and $\phi_E$
varies from 0 to 2. In this case the solution is simply that
$\phi_B=-2\phi_E$ as can be readily obtained from Equation
\eqref{eq:Phi_example} by setting $\Phi=\phi_A=0$.

\begin{figure}[htbp]
  \centering
  \SubFig{eq_phi_cm}{$\phi_A$, $\phi_B$}{\FigWidth}
  \SubFig{eq_x_cm}{$X$}{\FigWidth}
  \caption{Electrogenic system: Equilibria, order-one moieties. E is a
  chemostat, A \& B form a conserved moiety.
  % 
  (a) The potentials $\phi_A$ and $\phi_B$ when A and B form a
  conserved moiety with $x_A+x_B=2$ and the electrical potential $\phi_E$ varies.
  % 
  (b) The corresponding states $x_A$ and $x_B$. Note that $x_A+x_B=2$
  for all flows.}
\label{fig:cm}
\end{figure}
Consider secondly the case where only E is a chemostat. Then $\Ncd$
\eqref{eq:dX_cs} is given by $N$ with the third row
deleted:
\begin{align}
 \Ncd &=
  \begin{pmatrix}
    -1&1&0
  \end{pmatrix}^T
\end{align}
One choice of $G$ is
\begin{align}
  G &=
  \begin{pmatrix}
    1&1&0\\
    0&0&1
  \end{pmatrix}^T
\end{align}
As there are two non-zero entries in the first row of $G$, there is a
order-one conserved moiety and so the method of
\S~\ref{sec:conserved-moieties} must be used.
In this case
\begin{equation}
  \Gcm =
  \begin{pmatrix}
    1&1&0
  \end{pmatrix}
\end{equation}
Deleting the second element of the first row of $G$ gives the same
linear equation as for the two-chemostat case; the difference is that
$\phi_A$ must be chosen to satisfy the conserved moiety
$x_A+x_B=x_{AB}$. In this example, choose $x_{AB}=2$ and thus
$\phicm=\phi_N \ln 2$.  Figure \ref{fig:cm} shown the results of the procedure of
\S~\ref{sec:conserved-moieties} where $\phi_A$, $\phi_B$ plotted
against $\phi_E$ when the potential chemostat A is fixed at $\phi_A=0$
and $\phi_E$ varies from 0 to 2. 


\section{Steady-states}
\label{sec:steady-states}
In this paper, a system is said to be in \emph{steady-state} if efforts
are such that all system states are constant. Thus in the context of
stoichiometric models:
\begin{equation}\label{eq:ss}
  \dX = \Ncd \VV = 0
\end{equation}
where $\Ncd$ is the \emph{chemodynamic} stoichiometric matrix
\citep{GawCra16}.  In the special case of equilibrium
\eqref{eq:equilibrium}, $\VV = 0$ and Equation \eqref{eq:ss} is
automatically satisfied.

Pathway analysis determines the set  of \emph{non-zero} flows such
that Equation \eqref{eq:ss} is satisfied. In particular such analysis
determines the pathway matrix $K_p$ such that flows $\VV$ given by 
\begin{align}
  \VV &= K_p v_p
\end{align}
satisfy Equation \eqref{eq:ss} for any choice of $v_p$.


% But, in general, the flows $\VV$ are a \emph{nonlinear} function of
% potentials $\phi$.

%  which, following \S~\ref{sec:reaction-kinetics} can be
% \emph{approximately} written as a function of the reaction affinity
% $\Phi$ as 
% \begin{equation}\label{eq:F}
%   V = F(\Phi)
% \end{equation}
% where the function $F$ is the vector version of either the $\sinh$
% function of eqn.~\eqref{eq:MA_sinh_a} or the $\tanh$ function of
% eqn.~\eqref{eq:MA_tanh_a}.

% \begin{align}
%   \VV &= F_{fr}(\Phi_{fr}) \label{eq:F_fr}\\
%   \text{where } \Phi_{fr} &= 
%                             \begin{pmatrix}
%                               \Phi_f \\
%                               \hline
%                               \Phi_r 
%                             \end{pmatrix}
%   = -\NN_{fr}^T \phi\\
%   \text{and } \NN_{fr} &= 
%                           \begin{pmatrix}
%                             \NN_f \\
%                             \hline
%                             \NN_r 
%                           \end{pmatrix}
% \end{align}
% and $\Phi_{f}$ and $\Phi_{r}$ are the forward and reverse reaction
% affinities.    Hence solving Equations \eqref{eq:ss} and
% \eqref{eq:F_fr} for $\phi$ is a nonlinear problem. Some old, and a new,
% solutions to this problem are given in \S~\ref{sec:steady-states}.

% In some formulations, the  flows are a function of the net affinity:
% $\Phi = \Phi_{f}-\Phi_{r} = -N^T\phi$:
%   \begin{align}
%      \VV &= F(\Phi) \label{eq:F}
%   \end{align}




% \subsection{A (new?) positive-pathway approach}
% \label{sec:posit-pathw-appr}
\citet{GawCra17} present an energy-based approach to pathway
analysis. This forms the basis of the method discussed in this
section. In particular, the starting point for the method proposed
here is to apply positive-pathway analysis \citep{GawCra17} to
the system where the chemostats have been chosen to give exactly one
pathway described by $K_p$.

The method is as follows:
\begin{enumerate}
\item Choose a vector of pathway flows $v_p$ and deduce all
  the relevant steady-state flows from:
  \begin{equation}
    \VV = K_p v_p
  \end{equation}
\item Choose a starting value for the potentials $\phi$ -- for example
  that corresponding to equilibrium -- for all flows.
\item Compute the $\skappa$ -- based on $\phi$ -- for all flows.
\item \label{item:start} Deduce the reaction affinities by inverting  the approximation \eqref{eq:MA_sinh_a}
  \begin{equation}
    \Phi = 2\phi_N \asinh \frac{\VV}{2\kappabar}
  \end{equation}
\item Using the approach of \S~\ref{sec:equilibria}, deduce the
  free-to-vary potentials from Equation \eqref{eq:phi_cm_0}:
\begin{equation}
  \NN_v \phi_v = \Phi-\NN_{cs} \phi_{cs}
\end{equation}
% \item Use the method of \S~\ref{sec:conserved-moieties} to handle
%   conserved moieties where needed.
\item If $\phi$ has not converged, then return to step \ref{item:start}.
\end{enumerate}

Assuming convergence, this method yields the steady-state potentials
$\phi$ as well as the corresponding $\skappa$.

An alternative approach would use the $\tanh$ formulation
\eqref{eq:MA_tanh} and $\svv$ in place of the $\sinh$ formulation
\eqref{eq:MA_sinh} and $\skappa$.

% \subsection{Example: An extended Electrogenic Network}
% \label{sec:example:-an-extended}
% \begin{figure}[htbp]
%   \centering
%   \SubFig{ABCE_cbg}{Bond Graph}{0.45}
%   \SubFig{ABCE_phi_E}{$\phi_E$}{0.45}
%   \SubFig{ABCE_kappa}{$\dfrac{\skappa}{\skappaeq}$}{0.45}
%   \SubFig{ABCE_crite}{Convergence}{0.45}
%   \caption[An extended Electrogenic Network]{An Extended Electrogenic
%     Network.  
%     % 
%     (a) The system bond graph now has two reactions and an additional substance C.
%     % 
%     (b) The capacitor potential
%     plotted against reaction flow when $\phi_A=\SI{500}{mV}$ and
%     $\phi_B=\SI{100}{mV}$. 
%     % 
%     The reaction constants are again (see Figure
%     \ref{fig:RK-electrogenic}) chosen as 
%     % 
%     $\kappa_i = \frac{1}{\kappa_i^\star} = 8.4908e-06$.
%     % 
%     At $v=0$, $\phi_E = \frac{\phi_A - \phi_B}{2}=\SI{200}{mV}$; as
%     $v$ increases, $\phi_E$ decreases as the inverse hyperbolic sin:
%     $\asinh$. (c) $\skappa$ is \emph{not} constant. (d) Shows the
%     convergence of the method. }
% \end{figure}
% The example of \S~\ref{sec:exampl-basic-electr} is extended to a two-reaction case where the reaction
% \begin{equation}
%   A \reacu{r1} C \reacu{r2} B
% \end{equation}
% accumulates charge on the capacitor \BC{E}. The corresponding reaction
% affinities are:
% \begin{xalignat}{2}
%   \Phi^f_1 &= \phi_A - \phi_E & \Phi^f_2 &= \phi_C\\
%   \Phi^r_1 &=  \phi_C & \Phi^r_2 &= \phi_B + \phi_E\\
%   \Phi_1 &= \phi_A - \phi_C - \phi_E & \Phi_2 &= \phi_C - \phi_B - \phi_E\\
%   \Phi^{av}_1 &= \frac{\phi_A + \phi_C - \phi_E}{2} &
%   \Phi^{av}_2 &= \frac{\phi_B + \phi_C + \phi_E}{2}
% \end{xalignat}
% In the steady-state, both reaction flows are equal and given by:
% \begin{align}
%   v  &= 2 \kappa_1 \skappa_1 \sinh \frac{\phi_A - \phi_C -
%        \phi_E}{2\phi_N}\\
%      &= 2 \kappa_2 \skappa_2 \sinh \frac{\phi_C - \phi_B -
%        \phi_E}{2\phi_N}\\
% \text{where } \skappa_1 &= \exp \frac{\phi_A + \phi_B-
%                               \phi_E}{2}\\
% \text{and } \skappa_2 &= \exp \frac{\phi_B + \phi_C + \phi_E}{2}
% \end{align}
% In this case, $\skappa_1$ and $\skappa_2$ are \emph{not}
% constant and there is no direct solution of these two equations for
% $\phi_C$ and $\phi_E$. Instead, the iterative procedure of
% \S~\ref{sec:steady-states} is used. Figure \ref{subfig:ABCE_phi_E}
% shows how $\phi_E$ varies with $v$ after 50 iterations and Figure
% \ref{subfig:ABCE_kappa} shows how $\skappa$ 
% varies with $v_p$. Figure \ref{subfig:ABCE_crite} shows the convergence of
% the iterative approach.

% \subsection{Example: An Electrogenic Network with Conserved Moiety}
% \label{sec:example:-an-extended}
% \begin{figure}[htbp]
%   \centering
%   \SubFig{ABEFG_cbg}{Bond Graph}{0.45}
%   \SubFig{ABEFG_phi_E}{$\phi_E$}{0.45}
%   %%\SubFig{ABEFG_kappa}{$\dfrac{\skappa}{\skappaeq}$}{0.45}
%   \SubFig{ABEFG_crite}{Convergence}{0.45}
%   \SubFig{ABEFG_Xcm}{Conserved Moiety}{0.45}
%   \caption[Example: An Electrogenic Network with Conserved Moiety]
%   {Example: An Electrogenic Network with Conserved Moiety.
%     % 
%     (a) The system bond graph now has two reactions and two additional
%     substances F and G forming a conserved moiety. The initial states
%     are $x_F=x_G=1$.
%     % 
%     (b) The capacitor potential
%     plotted against reaction flow when $\phi_A=\SI{500}{mV}$,
%     $\phi_B=\SI{100}{mV}$ and $\phi_F=\phi_G=\SI{-100}{mV}$. 
%     % 
%     The first reaction constant is again (see Figure
%     \ref{fig:RK-electrogenic}) chosen as 
%     % 
%     $\kappa_1 = \frac{1}{\kappa_1^\star} = 8.4908e-06$.
%     % 
%     The second reaction constant is $\kappa_2 = \frac{200}{K_F+K_G}$.
%     % 
%     At $v=0$, $\phi_E = \frac{\phi_A - \phi_B}{2}=\SI{200}{mV}$; as
%     $v$ increases, $\phi_E$ decreases as the inverse hyperbolic sin:
%     $\asinh$. (c) $\skappa$ is \emph{not} constant. 
%     % 
%     (c) Shows the convergence of the method. 
%     % 
%     (d) Shows the two states $x_F$ and $x_G$ as reaction flow $v_p$
%     changes; $v_p=200$ is the maximum flow possible as $x_G<0$ for $v_p>200$.
% }\label{fig:ABEFG}
% \end{figure}
% The example of \S~\ref{sec:exampl-basic-electr} is extended to a
% two-reaction case where with a conserved moiety. The two reactions are:
% \begin{align}
%   A+G &\reacu{r1} B+F \\
%   F &\reacu{r2} G 
% \end{align}
% and charge is accumulated on the capacitor \BC{E}.

% Using the formulation of Equation \eqref{eq:K}, the second reaction
% flow can be written as:
% \begin{equation}\label{eq:v_2}
%   v_2 = \kappa_2 \lb K_Fx_F - K_G x_G\rb
% \end{equation}
% Using the conserved moiety: $x_F + x_G = x_{FG}$ where $x_{FG}$ is a
% constant, it follows from \eqref{eq:v_2} that:
% \begin{equation}
%   x_G = \frac{K_f}{K_F+K_G}x_{FG} - \frac{1}{\kappa_2\lb K_F+K_G\rb} v_2
% \end{equation}
% Thus $x_G$ decreases linearly with $v_2$. As $x_G\ge0$, the maximum
% possible flow $v_{max}$ is thus
% \begin{equation}
%   v_{max} = \kappa_2 K_F x_{FG}
% \end{equation}
% The results in Figure \ref{fig:ABEFG} are obtained for the special case
% that $\kappa_2 = \frac{200}{K_F+K_G}$ and $x_{FG}=2$ giving
% $v_{max}=200$.

% % \section{Example: A Simple Biomolecular Network}
% % \label{sec:exampl-simple-biom}
% % \begin{figure}[htbp]
% %   \centering
% %   \SubFig{Branch_cbg}{Bond Graph}{0.45}
% %   \SubFig{Branch_phi_P}{$\phi_P$}{0.45}
% %   \SubFig{Branch_phi_B}{$\phi_B$}{0.45}
% %   \SubFig{Branch_v3}{$v_3$}{0.45}
% %   \SubFig{Branch_Phi_av}{Normalised $\Phi_{av}$}{0.45}
% %   \caption{A Simple Biomolecular Network. $\kappa_1= 1,100,1000$, $\kappa_2=1$,
% % $\phi_A = 2$, $\phi_C=1$, $\phi_D=0.1$}
% % \end{figure}
% % Figure \ref{subfig:Branch_cbg} shown the bond graph of the following
% % simple biomolecular network:
% % \begin{align}
% %   A &\reacu{r1} B+P\\
% %   B &\reacu{r2} C+P\\
% %   P &\reacu{r3} D
% % \end{align}
% % This can be thought of as a simplified mitochondrial electron
% % transport chain where the two reactions $r1$ and $r2$ pump protons
% % into component \BC{P} and reaction $r3$ represents a leak of protons
% % to \BC{D}.

% % Assume initially that $v_3=0$.
% % At steady state, the potentials $\phi_A$, $\phi_B$, $\phi_c$ and
% % $\phi_P$ are related to the affinities $\Phi_1$ and $\Phi_2$ by:
% % \begin{align}
% %   \phi_A &= \phi_B + \phi_P + \Phi_1\\
% %   \phi_B &= \phi_C + \phi_P + \Phi_2
% % \end{align}
% % Solving these two simultaneous equations for $\phi_P$ and $\phi_B$ gives:
% % \begin{align}
% %   \phi_P &= \frac{\phi_A-\phi_C}{2} - \frac{\Phi_1+\Phi_2}{2}\\
% %   \phi_B &= \frac{\phi_A+\phi_C}{2} - \frac{\Phi_1-\Phi_2}{2}
% % \end{align}
% % Using the $\sinh$
% % approximation gives:
% % \begin{align}
% %   \Phi_1 &= 2 \phi_N \asinh \frac{v_p}{2\kappa_1}\\
% %   \Phi_2 &= 2 \phi_N \asinh \frac{v_p}{2\kappa_2}\\
% % \end{align}
% % where $v_p$ is the pathway flow.
% % In the particular case that $\kappa_1=\kappa_2=\kappa$,
% % \begin{align}
% %   \phi_P &= \frac{\phi_A-\phi_C}{2} - 2 \phi_N \asinh \frac{v_p}{2\kappa}\\
% %   \phi_B &= \frac{\phi_A+\phi_C}{2}
% % \end{align}
% % Figures \ref{subfig:Branch_phi_P}\&\ref{subfig:Branch_phi_B} show
% % $\phi_P$ and $\phi_B$ respectively for $-200 \le v_p \le 1000$. Note
% % the characteristic $\asinh$ shape of  Figure \ref{subfig:Branch_phi_P}.

\section{A Biomolecular Cycle}\label{sec:exampl-biom-cycle}
\begin{figure}[htbp]
  \centering
  \SubFig{Hill_diagram}{Schematic}{0.5}
  \SubFig{Hill_abg}{Bond Graph}{\FigWidth}
  %%\SubFig{HillE_abg}{Bond Graph: electrogenic}{0.35}
  \caption[Example: A Biomolecular Cycle]{Example: A Biomolecular
    Cycle.
    % 
    (a) Schematic corresponding to \citep[Figure 1.1]{Hil89}.
    % 
    (b) Bond graph  corresponding to (a).
    % 
    (c) Bond graph  with the electrogenic capacitor \BC{V_m} corresponding
    to a trans-membrane voltage.
}\label{fig:Cycle}
\end{figure}
This example is based on that of \citet[\S~5]{GawCra17} which in turn
corresponds to the seminal monograph, ``Free energy transduction and
biomolecular cycle kinetics'' of \citet{Hil89}. As discussed by
\citet[\S~5]{GawCra17}, the biomolecular cycle of Figure 1.1 of
\citet{Hil89} may be represented as the bond graph of Figure
\ref{subfig:Hill_abg}. 
% The corresponding chemical equations are:
% \begin{xalignat*}{2}
%    {M}_i + {E} &\reacu{em} {EM} &
%    {L}_i + {EM} &\reacu{lem} {LEM}\\
%    {LEM} &\reacu{lesm} {LE^\star M}&
%    {LE^\star M} &\reacu{esm} {L}_o + {E^\star M}\\
%    {E^\star M} &\reacu{es} {M}_o + {E^\star}&
%    {E^\star} &\reacu{e} {E}
%  \end{xalignat*}
\begin{figure}[htbp]
  \centering
  \SubFig{Hill_phi}{Pathway flow $v_p$ and $\phi_{Mi}\si{~mV}$}{\FigWidth}
  %%\SubFig{Hill_X}{Pathway flow $v_p$ and $x_{Mi}$}{0.3}
  \SubFig{Hill_X_E}{Conserved moiety}{0.45}
  \SubFig{Hill_rho}{$\rho_s$}{0.45}
  \caption[Example: A Biomolecular Cycle.]{Example: A
    Biomolecular Cycle.
    % 
    (a) The flow $v_P$ is plotted against the variable potential
    $\phi_{Mi}\si{~mV}$. The maximum flow $v_{max}$ is marked as a
    horizontal line.
    % 
    % (b) The flow $v_P$ is plotted against the variable state
    % $x_{Mi}$.
    % 
    (b) The six states corresponding to the conserved moiety are
    plotted against the pathway flow $v_P$. The maximum flow $v_{max}$
    is marked by a vertical line an corresponds to the flow where a
    state ($x_E$ in this case) becomes zero.
    % 
    (c) The values of $\rho_s$ \eqref{eq:MA_sinh} for the six reaction
    are, by definition, unity at zero flow (equilibrium), but are not
    unity for non-zero
    flow $v_p$.
}\label{fig:Hill}
\end{figure}
For the purposes of illustration the thermodynamic constants of the
ten species are $\phi^\std=0$ and $X^\std=1$ and the rate constants of
the six reactions are $\kappa=10$.
%
This system has a conserved moiety of order 5 involving the the
total amount of the six species ${E}$, ${EM}$, ${LEM}$, ${E^\star}$,
${E^\star M}$ and ${LE^\star M}$ and therefor the method of
\S~\ref{sec:conserved-moieties} is used.
%
For the purposes of illustration, the total amount is taken as
${E_{tot}=10}$.
%
Following \citep{GawCra17}, ${M}_o$, ${L}_i$, ${L}_o$ and ${M}_i$ are
chemostats and the amount of ${M}_o$ and ${L}_i$ is taken as unity,
the amount of ${L}_o$ as two and ${M}_i$ is variable.
%
With these parameters, as derived in the Supplementary Material of
\citep{GawCra17}, the maximum flow is given by
$v_{max}=\frac{100}{21} = 4.76$.

Figure \ref{subfig:Hill_phi} shows flow $v_P$ is plotted against the
variable potential $\phi_{Mi}\si{~mV}$. The maximum flow $v_{max}$ is
marked as a horizontal line.
%
Figure \ref{subfig:Hill_X_E} shows the amounts of the six species
involved in the conserved moiety; the sum remains 10 but one of the
amounts reaches zero at the maximum flow: this is  why there is a
maximum flow.
%
Figure \ref{subfig:Hill_rho} shows the values of $\rho_s$
\eqref{eq:MA_sinh} for the six reactions. By definition, $\rho_s=1$ at
zero flow (equilibrium), but is not unity for non-zero flow
$v_p$. Again, one of the $\rho_s$ reaches zero at the maximum flow.
% \begin{figure}[htbp]
%   \centering
%   \SubFig{HillE_phi}{$\phi_{Vm}\si{~mV}$}{\FigWidth}
%   %%\SubFig{HillE_X}{$x_{Vm}$}{\FigWidth}
%   \SubFig{HillE_X_E}{Conserved moiety}{\FigWidth}
%   \caption[Example: An Electrogenic Biomolecular Cycle]
%   {Example: An Electrogenic Biomolecular Cycle. Similar to Figure
%     \ref{fig:Hill} except that $x_{Mi}$ is constant and the membrane
%     potential $\phi_{Vm}\si{~mV}$ varies.
%   }\label{fig:HillE}
% \end{figure}

% Figure \ref{subfig:HillE_abg} is an illustrative electrogenic version
% of the bond graph of Figure \ref{subfig:Hill_abg} where with the
% electrogenic capacitor \BC{V_m} corresponds to a trans-membrane
% voltage. For this example, $x_{Mi}$ is constant and the membrane
% potential $\phi_{Vm}\si{~mV}$ varies. The results appear in Figure
% \ref{fig:HillE}.

\section{Mitochondrial Electron Transport Chain}
\label{sec:exampl-oxid-phosph}
% \begin{figure}[htbp]
%   \centering
%   \Fig{OxPhos_cbg}{1}
%   \caption{Oxidative phosphorylation. \textbf{mOx} and \textbf{mPhos}
%     are the modular versions of the Electron Transport Chain and the
%     ATPase phosphorylation reactions respectively.
%   }
% \end{figure}
\begin{figure}[htbp]
  \centering
  \Fig{Ox_cbg}{1}
  \caption{Mitochondrial Electron Transport Chain
    % 
    Following \citep{Gaw17a}, the three complexes CI, CIII and CIV
    are represented by the bond graph  modules \textbf{mCI},
    \textbf{mCIII} and \textbf{mCIV} (Figures 5, 6 and 7 of
    \citep{Gaw17a}) interconnected by energy bonds.
}
\label{fig:Ox}
\end{figure}

\begin{figure}[htbp]
  \centering
  \Fig{Ox_pmf}{\FigWidth}
  %%\SubFig{Ox_QFe}{Affinity Equalisation}{\FigWidth}
  \caption{Electron Transport Chain: PMF as flow varies.
    % 
    The firm line shows the results of steady-state analysis using the
    parameter values discussed in the text. The three points show the
    (normalised) experimetal data of \citet[Figure 2B (low \ch{Pi}
    data)]{BazBeaVin16} and the dashed line the simulation results
    that they present. As discussed in the text, a single parameter was
    adjusted to make the firm line fit the data.  }
  \label{fig:Ox_pmf}
\end{figure}
% \begin{figure}
%   \centering
%   \SubFig{OxPhos_pmf}{PMF $\psi$}{\FigWidth}
%   \SubFig{OxPhos_atp}{ATP potential}{\FigWidth}
%   \caption{Oxidative Phosphorylation}
% \end{figure}
As discovered by \citet{Mit93}, the key feature of
mitochondria is the \emph{chemiosmotic} energy transduction whereby a
chain of redox reactions pumps protons across the mitochondrial inner
membrane to generate the \emph{proton-motive force} (PMF). This PMF is
then used to power the synthesis of ATP -- the universal fuel of
living systems. The bond graph  model of chemiosmosis of Figure \ref{fig:Ox} is described in detail by
\citet{Gaw17a} and its equilibrium properties analysed.
%
This paper uses this model to illustrate the steady
state analysis of this biologically significant system and compares
the theoretical results to experimental results of
\citet{BazBeaVin16}.

\citet{BazBeaVin16} compare experimental results with a detailed ode
model of mitochondria which includes not only the Electron Transport
Chain but also the tricarboxylic acid (TCA) cycle and ATP
generation. The model is simulated to steady-state for a variety of
conditions and a number of model outputs are compared with
experimental data. One of these outputs is the membrane potential
$\psi$ which is compared with experimental data in \citep[Figure
2B]{BazBeaVin16} (Two sets of data are given: for low inorganic
phosphate (\ch{Pi}) and high \ch{Pi} -- the low \ch{Pi} data is used
here).  The values of intermembrane and matrix pH were used to compute
the corresponding proton-motive force (PMF).

 % ===== Discuss \citet{BazBeaVin16} here ====
% \begin{itemize}
% \item Includes TCA, ETC, ATPase
% \item Dynamic sim to steady state
% \item Compare varioius model outputs including $\psi$.
% \item This paper just looks at $\psi$ augmented to give pmf.
% \item 
% $\pH_x=7.2372$, $\pH_i = 7.2$, $\Phi_H = 2.2mV$, $\Delta_p = \psi + \Phi_H$
% \end{itemize}


The bond graph~\citep{Gaw17a} involves 10 chemical species: \ch{NADH},
\ch{NAD}, \ch{O2}, \ch{H2O}, \ch{Hx}, \ch{Hi}, \ch{Q}, \ch{QH2},
\ch{Fe3}, \ch{Fe2},
%
the transmembrane potential $\psi$ and the electrical charges used in
modelling the redox reactions. \ch{Q} and \ch{QH2} correspond to
oxidised and reduced ubiquinone respectively and \ch{Fe3} and \ch{Fe2}
correspond to oxidised and reduced cytochrome c respectively.
%
There are 9 reaction components\citep{Gaw17a}: the left \& right half
redox reactions and the proton pump corresponding to each of the three
complexes CI, CIII and CIV.
%

%
Following the analysis of \citet{GawCra17}, the positive pathway matrix is:
\begin{equation}
  K_p =
  \begin{pmatrix}
    2&2&4&2&4&4&4&1&4
  \end{pmatrix}^T
\end{equation}

%  The overall pathway equation is
% \begin{equation}
%   \ch{2 NADH + O2 + 18 Hx <> 2 NAD + 2 H2O + 16 Hi + 20 P}
% \end{equation}
% % \begin{equation}
%   \ch{2 NADH + O2 + 4 Hi <> 2 NAD + 2 H2O + 2 Hx + 20 PMF}
% \end{equation}

The corresponding reaction rates $\kappa_s$ of Equation
\eqref{eq:MA_sinh_0} are denoted $\kappa_{iL}$, $\kappa_{iR}$ and
$\kappa_{iP}$ coresponding to the left \& right half reactions and the
proton pump of the three complexes \citep{Gaw17a}. For the steady-state
analysis of Figure \ref{fig:Ox_pmf},
$\kappabar_{iL} = \kappa_{iR} = \infty$, $i=1\dots 3$ and
$\kappabar_{2P} = \kappabar_a$,
$\kappabar_{1P} = \kappabar_{3P}  =\infty$ where $\kappabar_a$ is
the sole adjustable parameter.  

As indicated in Figure \ref{fig:Ox_pmf}, when the sole adjustable
parameter $\kappabar_a=0.025$, this simple choice of adjustable
parameters, together with the Faraday-equivalent potentials listed in
\citep[Table 1]{Gaw17a}, gives a steady-state solution consitent with
the experimental data of \citet{BazBeaVin16}.

\section{Conclusion}
\label{sec:conclusion}
The mass-action reaction kinetics equation -- expressed in terms of
forward and reverse reaction affinity -- is rewritten in terms of
reaction affinity and summed reaction affinity. Combining this with
stoichiometric pathway analysis~\citep{GawCra17} leads to an algorithm
for computing steady state potentials corresponding to non-zero
steady-state reaction flows. This is illustrated using introductory
examples, a bond graph model of a biomolecular cycle, and a bond graph
model of the mitochondrial electron transport chain.

%%Future work:
Future work will extend the method to include non mass-action kinetics
such as Michaelis-Menten and to handle more than one pathway. It would
also be of interest to examine the relative merits of the $\sinh$ and
$\tanh$ and to analyse and enhance numerical properties of the
approach.
%

As noted in the Introduction, the approach in this paper is
specifically tailored to the special properties of biomolecular system
bond graphs. It would be interesting to see to what extent this
specific approach could be applied to the general bond graph approach
\citep{Bre84b,BidMarSca07,GonGal11}; similarly it would be interesting
to see how derivative causality appled to biomolecular system bond
graphs \citep{Gaw17} could be combined with the general approach
\citep{Bre84b,BidMarSca07,GonGal11} to deduce structural properties
from the bond graph.

The illustrative example of \S~\ref{sec:exampl-oxid-phosph} just
matches the proton-motive force (PMF) to the
data~\citep{BazBeaVin16}. It would be interesting to use the other
free parameters to match the other types of data analysed by
\citet{BazBeaVin16}.


% \begin{enumerate}
% \item Numerical analysis of the iterative algorithm and application of
%   numerical methods to improve convergence properties.
% % \item Use the method to fit parameters to the bond graph  model of the
% %   ETC to fit \emph{all} the data given in Figure 2
% %   \citet{BazBeaVin16}; at the moment it just fits Figure 2B: the PMF
% %   $\psi$.
% \item Compare and contrast the $\sinh$~\eqref{eq:MA_sinh_a} and
%   $\tanh$~\eqref{eq:MA_tanh_a} formulations.
% \item Extend to Michaelis-Menten reaction kinetics. 
% \item Examine examples with more than one independent pathway.
% \item Apply the approach to the bond graph  model of the mitochondrial
%   electron transport chain derived by \citep{Gaw17a} and fit to the
%   experimental results of \citep{BazBeaVin16}.
% \end{enumerate}

\section{Acknowledgements}
Peter Gawthrop would like to thank the Melbourne School of Engineering
for its support via a Professorial Fellowship, Edmund Crampin for
help, advice and encouragement and Michael Pan for comments on an
earlier draft.
%
The constructive input of the reviewers is greatfully acknowledged.

%%\bibliography{common}
\begin{thebibliography}{29}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{}
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{doi: #1}\else
  \providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi

\bibitem[Beard(2011)]{Bea11}
Daniel~A. Beard.
\newblock Simulation of cellular biochemical system kinetics.
\newblock \emph{Wiley Interdisciplinary Reviews: Systems Biology and Medicine},
  3\penalty0 (2):\penalty0 136--146, 2011.
\newblock ISSN 1939-005X.
\newblock \doi{10.1002/wsbm.116}.

\bibitem[Beard(2012)]{Bea12}
Daniel~A. Beard.
\newblock \emph{{Biosimulation: Simulation of Living Systems}}.
\newblock Cambridge University Press, Cambridge, UK., 2012.
\newblock ISBN 978-0-521-76823-8.

\bibitem[Beard and Qian(2010)]{BeaQia10}
Daniel~A Beard and Hong Qian.
\newblock \emph{Chemical biophysics: quantitative analysis of cellular
  systems}.
\newblock Cambridge University Press, 2010.

\bibitem[Atkins and de~Paula(2011)]{AtkPau11}
Peter Atkins and Julio de~Paula.
\newblock \emph{{Physical Chemistry for the Life Sciences}}.
\newblock Oxford University Press, 2nd edition, 2011.

\bibitem[Oster et~al.(1971)Oster, Perelson, and Katchalsky]{OstPerKat71}
George Oster, Alan Perelson, and Aharon Katchalsky.
\newblock Network thermodynamics.
\newblock \emph{Nature}, 234:\penalty0 393--399, December 1971.
\newblock \doi{10.1038/234393a0}.

\bibitem[Oster et~al.(1973)Oster, Perelson, and Katchalsky]{OstPerKat73}
George~F. Oster, Alan~S. Perelson, and Aharon Katchalsky.
\newblock Network thermodynamics: dynamic modelling of biophysical systems.
\newblock \emph{Quarterly Reviews of Biophysics}, 6\penalty0 (01):\penalty0
  1--134, 1973.
\newblock \doi{10.1017/S0033583500000081}.

\bibitem[Greifeneder and Cellier(2012)]{GreCel12}
J.~Greifeneder and F.E. Cellier.
\newblock Modeling chemical reactions using bond graphs.
\newblock In \emph{Proceedings ICBGM12, 10th SCS Intl. Conf. on Bond Graph
  Modeling and Simulation, Genoa, Italy}, pages 110--121, San Diego, CA, USA,
  July 2012. The Society for Modeling and Simulation International.

\bibitem[Gawthrop and Crampin(2014)]{GawCra14}
Peter~J. Gawthrop and Edmund~J. Crampin.
\newblock Energy-based analysis of biochemical cycles using bond graphs.
\newblock \emph{Proceedings of the Royal Society A: Mathematical, Physical and
  Engineering Science}, 470\penalty0 (2171):\penalty0 1--25, 2014.
\newblock \doi{10.1098/rspa.2014.0459}.
\newblock Available at {arXiv:1406.2447}.

\bibitem[Gawthrop and Crampin(2016)]{GawCra16}
P.~J. Gawthrop and E.~J. Crampin.
\newblock Modular bond-graph modelling and analysis of biomolecular systems.
\newblock \emph{IET Systems Biology}, 10\penalty0 (5):\penalty0 187--201,
  October 2016.
\newblock ISSN 1751-8849.
\newblock \doi{10.1049/iet-syb.2015.0083}.
\newblock Available at {arXiv:1511.06482}.

\bibitem[Gawthrop and Crampin(2017)]{GawCra17}
Peter~J. Gawthrop and Edmund~J. Crampin.
\newblock Energy-based analysis of biomolecular pathways.
\newblock \emph{Proceedings of the Royal Society of London A: Mathematical,
  Physical and Engineering Sciences}, 473\penalty0 (2202), 2017.
\newblock ISSN 1364-5021.
\newblock \doi{10.1098/rspa.2016.0825}.
\newblock Available at {arXiv:1611.02332}.

\bibitem[Gawthrop et~al.(2017)Gawthrop, Siekmann, Kameneva, Saha, Ibbotson, and
  Crampin]{GawSieKam17}
P.~J. Gawthrop, I.~Siekmann, T.~Kameneva, S.~Saha, M.~R. Ibbotson, and E.~J.
  Crampin.
\newblock Bond graph modelling of chemoelectrical energy transduction.
\newblock \emph{IET Systems Biology}, 11\penalty0 (5):\penalty0 127--138, 2017.
\newblock ISSN 1751-8849.
\newblock \doi{10.1049/iet-syb.2017.0006}.
\newblock Available at {arXiv:1512.00956}.

\bibitem[Gawthrop(2017{\natexlab{a}})]{Gaw17a}
P.~J. Gawthrop.
\newblock Bond graph modeling of chemiosmotic biomolecular energy transduction.
\newblock \emph{IEEE Transactions on NanoBioscience}, 16\penalty0 (3):\penalty0
  177--188, April 2017{\natexlab{a}}.
\newblock ISSN 1536-1241.
\newblock \doi{10.1109/TNB.2017.2674683}.
\newblock Available at {arXiv:1611.04264}.

\bibitem[Borutzky(2010)]{Bor10}
Wolfgang Borutzky.
\newblock \emph{Bond graph methodology: development and analysis of
  multidisciplinary dynamic system models}.
\newblock Springer, 2010.
\newblock ISBN 978-1-84882-881-0.
\newblock \doi{10.1007/978-1-84882-882-7}.

\bibitem[Karnopp et~al.(2012)Karnopp, Margolis, and Rosenberg]{KarMarRos12}
Dean~C Karnopp, Donald~L Margolis, and Ronald~C Rosenberg.
\newblock \emph{System Dynamics: Modeling, Simulation, and Control of
  Mechatronic Systems}.
\newblock John Wiley \& Sons, 5th edition, 2012.
\newblock ISBN 978-0470889084.

\bibitem[Gawthrop and Bevan(2007)]{GawBev07}
Peter~J Gawthrop and Geraint~P Bevan.
\newblock Bond-graph modeling: A tutorial introduction for control engineers.
\newblock \emph{IEEE Control Systems Magazine}, 27\penalty0 (2):\penalty0
  24--45, April 2007.
\newblock \doi{10.1109/MCS.2007.338279}.

\bibitem[Park et~al.(2016)Park, Rubin, Xu, Amador-Noguez, Fan, Shlomi, and
  Rabinowitz]{ParRubXu16}
Junyoung~O. Park, Sara~A. Rubin, Yi-Fan Xu, Daniel Amador-Noguez, Jing Fan,
  Tomer Shlomi, and Joshua~D. Rabinowitz.
\newblock Metabolite concentrations, fluxes and free energies imply efficient
  enzyme usage.
\newblock \emph{Nat Chem Biol}, 12\penalty0 (7):\penalty0 482--489, Jul 2016.
\newblock ISSN 1552-4450.
\newblock \doi{10.1038/nchembio.2077}.

\bibitem[Orth et~al.(2010)Orth, Thiele, and Palsson]{OrtThiPal10}
Jeffrey~D. Orth, Ines Thiele, and Bernhard~O. Palsson.
\newblock What is flux balance analysis?
\newblock \emph{Nat Biotech}, 28:\penalty0 245--248, March 2010.
\newblock ISSN 1087-0156.
\newblock \doi{10.1038/nbt.1614}.

\bibitem[King and Altman(1956)]{KinAlt56}
Edward~L. King and Carl Altman.
\newblock A schematic method of deriving the rate laws for enzyme-catalyzed
  reactions.
\newblock \emph{The Journal of Physical Chemistry}, 60\penalty0 (10):\penalty0
  1375--1378, 1956.
\newblock \doi{10.1021/j150544a010}.

\bibitem[Hill(1989)]{Hil89}
Terrell~L Hill.
\newblock \emph{Free energy transduction and biochemical cycle kinetics}.
\newblock Springer-Verlag, New York, 1989.

\bibitem[Bazil et~al.(2016)Bazil, Beard, and Vinnakota]{BazBeaVin16}
Jason~N. Bazil, Daniel~A. Beard, and Kalyan~C. Vinnakota.
\newblock Catalytic coupling of oxidative phosphorylation, {ATP} demand, and
  reactive oxygen species generation.
\newblock \emph{Biophysical Journal}, 110\penalty0 (4):\penalty0 962 -- 971,
  2016.
\newblock ISSN 0006-3495.
\newblock \doi{10.1016/j.bpj.2015.09.036}.

\bibitem[Henry et~al.(2007)Henry, Broadbelt, and Hatzimanikatis]{HenBroHat07}
Christopher~S. Henry, Linda~J. Broadbelt, and Vassily Hatzimanikatis.
\newblock Thermodynamics-based metabolic flux analysis.
\newblock \emph{Biophysical Journal}, 92\penalty0 (5):\penalty0 1792 -- 1805,
  2007.
\newblock ISSN 0006-3495.
\newblock \doi{10.1529/biophysj.106.093138}.

\bibitem[Qi et~al.(2009)Qi, Dash, Han, and Beard]{QiDasHan09}
Feng Qi, Ranjan~K. Dash, Yu~Han, and Daniel~A. Beard.
\newblock Generating rate equations for complex enzyme systems by a
  computer-assisted systematic method.
\newblock \emph{BMC Bioinformatics}, 10\penalty0 (1):\penalty0 1--9, 2009.
\newblock ISSN 1471-2105.
\newblock \doi{10.1186/1471-2105-10-238}.

\bibitem[Breedveld(1984)]{Bre84b}
P.C. Breedveld.
\newblock A bond graph algorithm to determine the equilibrium state of a
  system.
\newblock \emph{Journal of the Franklin Institute}, 318\penalty0 (2):\penalty0
  71 -- 75, 1984.
\newblock ISSN 0016-0032.
\newblock \doi{10.1016/S0016-0032(84)90006-1}.

\bibitem[Bideaux et~al.(2006)Bideaux, Marquis-Favre, and Scavarda]{BidMarSca07}
Eric Bideaux, Wilfrid Marquis-Favre, and Serge Scavarda.
\newblock Equilibrium set investigation using bicausality.
\newblock \emph{Mathematical and Computer Modelling of Dynamical Systems},
  12\penalty0 (2-3):\penalty0 127--140, 2006.
\newblock \doi{10.1080/13873950500069003}.

\bibitem[Gonz�É�lez and Galindo(2011)]{GonGal11}
G~Gonz�É�lez and R~Galindo.
\newblock Steady state determination using bond graphs for systems with
  singular state matrix.
\newblock \emph{Proceedings of the Institution of Mechanical Engineers, Part I:
  Journal of Systems and Control Engineering}, 225\penalty0 (7):\penalty0
  887--901, 2011.
\newblock \doi{10.1177/2041304110394552}.

\bibitem[Gawthrop(2017{\natexlab{b}})]{Gaw17}
Peter~J. Gawthrop.
\newblock Bond-graph modelling and causal analysis of biomolecular systems.
\newblock In Wolfgang Borutzky, editor, \emph{Bond Graphs for Modelling,
  Control and Fault Diagnosis of Engineering Systems}, pages 587--623. Springer
  International Publishing, Berlin, 2017{\natexlab{b}}.
\newblock ISBN 978-3-319-47434-2.
\newblock \doi{10.1007/978-3-319-47434-2_16}.

\bibitem[Polettini and Esposito(2014)]{PolEsp14}
Matteo Polettini and Massimiliano Esposito.
\newblock {Irreversible thermodynamics of open chemical networks. I. Emergent
  cycles and broken conservation laws}.
\newblock \emph{The Journal of Chemical Physics}, 141\penalty0 (2):\penalty0
  024117, 2014.
\newblock \doi{10.1063/1.4886396}.

\bibitem[Eaton et~al.(2015)Eaton, Bateman, Hauberg, and Wehbring]{EatBatHau15}
John~W. Eaton, David Bateman, S{\o}ren Hauberg, and Rik Wehbring.
\newblock \emph{{GNU Octave} version 4.0.0 manual: a high-level interactive
  language for numerical computations}, 2015.
\newblock URL \url{http://www.gnu.org/software/octave/doc/interpreter}.

\bibitem[Mitchell(1993)]{Mit93}
Peter Mitchell.
\newblock {David Keilin���s Respiratory Chain Concept and its Chemiosmotic
  Consequences}.
\newblock In Tore Fr{\"a}ngsmyr and Sture Fors{\'e}n, editors, \emph{Nobel
  Lectures in Chemistry, 1971-1980}. World Scientific, Singapore, 1993.
\newblock ISBN 981-02-0786-7.

\end{thebibliography}

\end{document}

%%% Local Variables:
%%% mode: latex
%%% TeX-master: t
%%% End:

1706.10278.extracted_bib
[Beard(2011)]{Bea11 Daniel~A. Beard. \newblock Simulation of cellular biochemical system kinetics. \newblock Wiley Interdisciplinary Reviews: Systems Biology and Medicine, 3\penalty0 (2):\penalty0 136--146, 2011. \newblock ISSN 1939-005X. \newblock \doi{10.1002/wsbm.116.
1706.10278.extracted_bib
[Beard(2012)]{Bea12 Daniel~A. Beard. \newblock {Biosimulation: Simulation of Living Systems. \newblock Cambridge University Press, Cambridge, UK., 2012. \newblock ISBN 978-0-521-76823-8.
1706.10278.extracted_bib
[Beard and Qian(2010)]{BeaQia10 Daniel~A Beard and Hong Qian. \newblock Chemical biophysics: quantitative analysis of cellular systems. \newblock Cambridge University Press, 2010.
1706.10278.extracted_bib
[Atkins and de~Paula(2011)]{AtkPau11 Peter Atkins and Julio de~Paula. \newblock {Physical Chemistry for the Life Sciences. \newblock Oxford University Press, 2nd edition, 2011.
1706.10278.extracted_bib
[Oster et~al.(1971)Oster, Perelson, and Katchalsky]{OstPerKat71 George Oster, Alan Perelson, and Aharon Katchalsky. \newblock Network thermodynamics. \newblock Nature, 234:\penalty0 393--399, December 1971. \newblock \doi{10.1038/234393a0.
1706.10278.extracted_bib
[Oster et~al.(1973)Oster, Perelson, and Katchalsky]{OstPerKat73 George~F. Oster, Alan~S. Perelson, and Aharon Katchalsky. \newblock Network thermodynamics: dynamic modelling of biophysical systems. \newblock Quarterly Reviews of Biophysics, 6\penalty0 (01):\penalty0 1--134, 1973. \newblock \doi{10.1017/S0033583500000081.
1706.10278.extracted_bib
[Greifeneder and Cellier(2012)]{GreCel12 J.~Greifeneder and F.E. Cellier. \newblock Modeling chemical reactions using bond graphs. \newblock In Proceedings ICBGM12, 10th SCS Intl. Conf. on Bond Graph Modeling and Simulation, Genoa, Italy, pages 110--121, San Diego, CA, USA, July 2012. The Society for Modeling and Simulation International.
1706.10278.extracted_bib
[Gawthrop and Crampin(2014)]{GawCra14 Peter~J. Gawthrop and Edmund~J. Crampin. \newblock Energy-based analysis of biochemical cycles using bond graphs. \newblock Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, 470\penalty0 (2171):\penalty0 1--25, 2014. \newblock \doi{10.1098/rspa.2014.0459. \newblock Available at {arXiv:1406.2447.
1706.10278.extracted_bib
[Gawthrop and Crampin(2016)]{GawCra16 P.~J. Gawthrop and E.~J. Crampin. \newblock Modular bond-graph modelling and analysis of biomolecular systems. \newblock IET Systems Biology, 10\penalty0 (5):\penalty0 187--201, October 2016. \newblock ISSN 1751-8849. \newblock \doi{10.1049/iet-syb.2015.0083. \newblock Available at {arXiv:1511.06482.
1706.10278.extracted_bib
[Gawthrop and Crampin(2017)]{GawCra17 Peter~J. Gawthrop and Edmund~J. Crampin. \newblock Energy-based analysis of biomolecular pathways. \newblock Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 473\penalty0 (2202), 2017. \newblock ISSN 1364-5021. \newblock \doi{10.1098/rspa.2016.0825. \newblock Available at {arXiv:1611.02332.
1706.10278.extracted_bib
[Gawthrop et~al.(2017)Gawthrop, Siekmann, Kameneva, Saha, Ibbotson, and Crampin]{GawSieKam17 P.~J. Gawthrop, I.~Siekmann, T.~Kameneva, S.~Saha, M.~R. Ibbotson, and E.~J. Crampin. \newblock Bond graph modelling of chemoelectrical energy transduction. \newblock IET Systems Biology, 11\penalty0 (5):\penalty0 127--138, 2017. \newblock ISSN 1751-8849. \newblock \doi{10.1049/iet-syb.2017.0006. \newblock Available at {arXiv:1512.00956.
1706.10278.extracted_bib
[Gawthrop(2017{\natexlab{a)]{Gaw17a P.~J. Gawthrop. \newblock Bond graph modeling of chemiosmotic biomolecular energy transduction. \newblock IEEE Transactions on NanoBioscience, 16\penalty0 (3):\penalty0 177--188, April 2017{\natexlab{a. \newblock ISSN 1536-1241. \newblock \doi{10.1109/TNB.2017.2674683. \newblock Available at {arXiv:1611.04264.
1706.10278.extracted_bib
[Borutzky(2010)]{Bor10 Wolfgang Borutzky. \newblock Bond graph methodology: development and analysis of multidisciplinary dynamic system models. \newblock Springer, 2010. \newblock ISBN 978-1-84882-881-0. \newblock \doi{10.1007/978-1-84882-882-7.
1706.10278.extracted_bib
[Karnopp et~al.(2012)Karnopp, Margolis, and Rosenberg]{KarMarRos12 Dean~C Karnopp, Donald~L Margolis, and Ronald~C Rosenberg. \newblock System Dynamics: Modeling, Simulation, and Control of Mechatronic Systems. \newblock John Wiley \& Sons, 5th edition, 2012. \newblock ISBN 978-0470889084.
1706.10278.extracted_bib
[Gawthrop and Bevan(2007)]{GawBev07 Peter~J Gawthrop and Geraint~P Bevan. \newblock Bond-graph modeling: A tutorial introduction for control engineers. \newblock IEEE Control Systems Magazine, 27\penalty0 (2):\penalty0 24--45, April 2007. \newblock \doi{10.1109/MCS.2007.338279.
1706.10278.extracted_bib
[Park et~al.(2016)Park, Rubin, Xu, Amador-Noguez, Fan, Shlomi, and Rabinowitz]{ParRubXu16 Junyoung~O. Park, Sara~A. Rubin, Yi-Fan Xu, Daniel Amador-Noguez, Jing Fan, Tomer Shlomi, and Joshua~D. Rabinowitz. \newblock Metabolite concentrations, fluxes and free energies imply efficient enzyme usage. \newblock Nat Chem Biol, 12\penalty0 (7):\penalty0 482--489, Jul 2016. \newblock ISSN 1552-4450. \newblock \doi{10.1038/nchembio.2077.
1706.10278.extracted_bib
[Orth et~al.(2010)Orth, Thiele, and Palsson]{OrtThiPal10 Jeffrey~D. Orth, Ines Thiele, and Bernhard~O. Palsson. \newblock What is flux balance analysis? \newblock Nat Biotech, 28:\penalty0 245--248, March 2010. \newblock ISSN 1087-0156. \newblock \doi{10.1038/nbt.1614.
1706.10278.extracted_bib
[King and Altman(1956)]{KinAlt56 Edward~L. King and Carl Altman. \newblock A schematic method of deriving the rate laws for enzyme-catalyzed reactions. \newblock The Journal of Physical Chemistry, 60\penalty0 (10):\penalty0 1375--1378, 1956. \newblock \doi{10.1021/j150544a010.
1706.10278.extracted_bib
[Hill(1989)]{Hil89 Terrell~L Hill. \newblock Free energy transduction and biochemical cycle kinetics. \newblock Springer-Verlag, New York, 1989.
1706.10278.extracted_bib
[Bazil et~al.(2016)Bazil, Beard, and Vinnakota]{BazBeaVin16 Jason~N. Bazil, Daniel~A. Beard, and Kalyan~C. Vinnakota. \newblock Catalytic coupling of oxidative phosphorylation, {ATP demand, and reactive oxygen species generation. \newblock Biophysical Journal, 110\penalty0 (4):\penalty0 962 -- 971, 2016. \newblock ISSN 0006-3495. \newblock \doi{10.1016/j.bpj.2015.09.036.
1706.10278.extracted_bib
[Henry et~al.(2007)Henry, Broadbelt, and Hatzimanikatis]{HenBroHat07 Christopher~S. Henry, Linda~J. Broadbelt, and Vassily Hatzimanikatis. \newblock Thermodynamics-based metabolic flux analysis. \newblock Biophysical Journal, 92\penalty0 (5):\penalty0 1792 -- 1805, 2007. \newblock ISSN 0006-3495. \newblock \doi{10.1529/biophysj.106.093138.
1706.10278.extracted_bib
[Qi et~al.(2009)Qi, Dash, Han, and Beard]{QiDasHan09 Feng Qi, Ranjan~K. Dash, Yu~Han, and Daniel~A. Beard. \newblock Generating rate equations for complex enzyme systems by a computer-assisted systematic method. \newblock BMC Bioinformatics, 10\penalty0 (1):\penalty0 1--9, 2009. \newblock ISSN 1471-2105. \newblock \doi{10.1186/1471-2105-10-238.
1706.10278.extracted_bib
[Breedveld(1984)]{Bre84b P.C. Breedveld. \newblock A bond graph algorithm to determine the equilibrium state of a system. \newblock Journal of the Franklin Institute, 318\penalty0 (2):\penalty0 71 -- 75, 1984. \newblock ISSN 0016-0032. \newblock \doi{10.1016/S0016-0032(84)90006-1.
1706.10278.extracted_bib
[Bideaux et~al.(2006)Bideaux, Marquis-Favre, and Scavarda]{BidMarSca07 Eric Bideaux, Wilfrid Marquis-Favre, and Serge Scavarda. \newblock Equilibrium set investigation using bicausality. \newblock Mathematical and Computer Modelling of Dynamical Systems, 12\penalty0 (2-3):\penalty0 127--140, 2006. \newblock \doi{10.1080/13873950500069003.
1706.10278.extracted_bib
[Gonz�É�lez and Galindo(2011)]{GonGal11 G~Gonz�É�lez and R~Galindo. \newblock Steady state determination using bond graphs for systems with singular state matrix. \newblock Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, 225\penalty0 (7):\penalty0 887--901, 2011. \newblock \doi{10.1177/2041304110394552.
1706.10278.extracted_bib
[Gawthrop(2017{\natexlab{b)]{Gaw17 Peter~J. Gawthrop. \newblock Bond-graph modelling and causal analysis of biomolecular systems. \newblock In Wolfgang Borutzky, editor, Bond Graphs for Modelling, Control and Fault Diagnosis of Engineering Systems, pages 587--623. Springer International Publishing, Berlin, 2017{\natexlab{b. \newblock ISBN 978-3-319-47434-2. \newblock \doi{10.1007/978-3-319-47434-2_16.
1706.10278.extracted_bib
[Polettini and Esposito(2014)]{PolEsp14 Matteo Polettini and Massimiliano Esposito. \newblock {Irreversible thermodynamics of open chemical networks. I. Emergent cycles and broken conservation laws. \newblock The Journal of Chemical Physics, 141\penalty0 (2):\penalty0 024117, 2014. \newblock \doi{10.1063/1.4886396.
1706.10278.extracted_bib
[Eaton et~al.(2015)Eaton, Bateman, Hauberg, and Wehbring]{EatBatHau15 John~W. Eaton, David Bateman, S{\oren Hauberg, and Rik Wehbring. \newblock {GNU Octave version 4.0.0 manual: a high-level interactive language for numerical computations, 2015. \newblock URL \url{http://www.gnu.org/software/octave/doc/interpreter.
1706.10278.extracted_bib
[Mitchell(1993)]{Mit93 Peter Mitchell. \newblock {David Keilin���s Respiratory Chain Concept and its Chemiosmotic Consequences. \newblock In Tore Fr{\"angsmyr and Sture Fors{\'en, editors, Nobel Lectures in Chemistry, 1971-1980. World Scientific, Singapore, 1993. \newblock ISBN 981-02-0786-7.