Physically-Plausible Modelling of Biomolecular Systems: A Simplified, Energy-Based Model of the Mitochondrial Electron Transport Chain

TL;DR

This paper introduces a simplified, thermodynamically consistent model of the mitochondrial electron transport chain that captures essential features without full mechanistic detail, aiding large-scale cellular modeling.

Contribution

The paper presents a novel approach for creating simplified, physically plausible models of complex biochemical systems suitable for integration into comprehensive cellular models.

Findings

The simplified model replicates key behaviors of the full electron transport chain.

Thermodynamic consistency is maintained in the simplified model.

The approach facilitates large-scale cellular biochemistry modeling.

Abstract

Systems biology and whole-cell modelling are demanding increasingly comprehensive mathematical models of cellular biochemistry. These models require the development of simplified models of specific processes which capture essential biophysical features but without unnecessarily complexity. Recently there has been renewed interest in thermodynamically-based modelling of cellular processes. Here we present an approach to developing of simplified yet thermodynamically consistent (hence physically plausible) models which can readily be incorporated into large scale biochemical descriptions but which do not require full mechanistic detail of the underlying processes. We illustrate the approach through development of a simplified, physically plausible model of the mitochondrial electron transport chain and show that the simplified model behaves like the full system.

Equations

$$\mathrm{ch}A<>[r]2B$$

$${\mu }_A$$

$${\mu }_B$$

$$f=\kappa \left(\exp \frac{A^f}{RT}-\exp \frac{A^r}{RT}\right)$$

$$A^f$$

$$\text{and }{A}^r$$

$$f=\kappa \left(K_A{x}_A-K_B^2{x}_B^2\right)$$

$$f_x$$

$$N={\left(\begin{array}{cc}-1& 2\end{array}\right)}^T$$

$$f=\kappa \left(\sqrt{K_A{x}_A}-K_B{x}_B\right)$$

$$N_{\alpha }={\left(\begin{array}{cc}-\frac{1}{2}& 1\end{array}\right)}^T$$

$$f$$

$$\text{where }{\kappa }_{\alpha }(x_A,x_B)$$

$$f$$

$$\varphi$$

$$f$$

$${\varphi }_A$$

$${\varphi }_B$$

$$\text{where }{V}_N$$

$$f$$

$$\mathrm{ch}A+2B<>[r]C+2D$$

$$\mathrm{ch}2NADH+O2+22H_x+<>2NAD++2H2O+20H_i+$$

Full text

1 Introduction

In mathematical biology, and more widely, the relative merits of simple ‘toy’ models, which represent some key aspects of the system but not full mechanistic detail, and comprehensive mechanistically detailed representations have long been debated. Simple ‘toy’ models allow rigorous mathematical analysis and are generally easy to simulate, but are difficult to relate to the full system and measurements thereof. Full mechanistically detailed models on the other hand provide a straight-forward mapping to the real system, but are challenging to parameterize and analyse, and may require significant computational overhead to simulate.

Simple models of complex biochemical processes can elucidate basic behaviour and biologically significant trade-offs (Scott et al., 2014; Weiße et al., 2015) and as such can be used as an aid to synthetic biology (Darlington et al., 2018). Furthermore, models of individual processes may be used as part of a model of an overall system as, for example, in the Physiome Project (Crampin et al., 2004; Hunter, 2016), or in whole-cell modelling Karr et al. (2012); Macklin et al. (2014); this requires models to be modular and reusable (Neal et al., 2014; Nickerson et al., 2016).

Recently there has been renewed interest in thermodynamically-based mechanistic modelling of cellular processes (Mason and Covert, 2019; Pan et al., 2019; Gawthrop et al., 2017; Klipp et al., 2016; Beard and Qian, 2010). A modular approach to energy-based modelling has been developed in the context of biomolecular systems (Gawthrop et al., 2015; Gawthrop and Crampin, 2016). This raises the question as to whether it is possible to develop energy-based models that are nevertheless simple.

Like engineering systems, living systems are subject to the laws of physics in general and the laws of thermodynamics in particular. This fact gives the opportunity of applying engineering approaches to the modelling, analysis and understanding of living systems. The bond graph method of Paynter (1961) is one such well-established engineering approach (Cellier, 1991; Gawthrop and Smith, 1996; Gawthrop and Bevan, 2007; Borutzky, 2010; Karnopp et al., 2012) which has been extended to include biomolecular systems (Oster et al., 1971, 1973; Gawthrop and Crampin, 2014, 2017; Gawthrop et al., 2017; Gawthrop and Crampin, 2018a, b; Pan et al., 2018, 2019).

When developing simplified models of biomolecular systems where energy transduction is important, it is essential that models be physically-plausible. A physically-plausible model of a physical system has two attributes: it is itself a model of a physical system (i.e. it does not contravene the laws of physics); and it shares key behaviours with the actual physical system (Gawthrop, 2003). Such an approach will, however, only be of use if there are complex physical systems which can indeed be represented by a simpler physical model. This paper shows that this is indeed the case. In particular we demonstrate that it is possible to develop a simplified model of the mitochondrial electron transport chain that is thermodynamically consistent, but which doesn’t represent full mechanistic detail, and show that it behaves like the full model.

Mitochondria make use of reduction-oxidation (redox) reactions in which the transfer of electrons is used to provide the power driving many living systems. As discovered by Mitchell (1961, 1976, 1993, 2011), the key feature of mitochondria is the chemiosmotic energy transduction whereby a chain of redox reactions pumps protons across the mitochondrial inner membrane to generate an electrochemical gradient known as the proton-motive force (PMF). The PMF is then used to power the synthesis of ATP – the universal fuel of living systems. Due to this central role in living systems, mathematical modelling of the key components of mitochondria is thus an important challenge to systems biology. Because mitochondria transduce energy, an energy-based modelling method is desirable, and Beard and colleagues have developed the most comprehensive such models to date (Beard, 2005; Wu et al., 2007; Beard and Qian, 2010; Beard, 2012; Bazil et al., 2016). A bond graph model of mitochondrial oxidative phosphorylation has been given by Gawthrop (2017). This model is based on modelling the redox reactions associated with complexes \chCI, \chCIII and \chCIV of the mitochondrial electron transport chain.

Below we briefly outline the bond graph approach to modelling energy flows in biochemical reactions, in particular describing the Faraday-equivalent potential approach to modelling electrochemical phenomena, and we describe a modified mass action kinetics approach which will be central to development of a simplified thermodynamic modelling approach. A set of Python based tools has been developed to assist the development and analysis of bond graph models and these tools are briefly outlined.

We then discuss the Mitochondrial Electron Transport Chain (ETC) as an example of a complex biomolecular system which can be successfully modelled by a simple physically-plausible model, and use data from Bazil et al. (2016) to derive parameters of the physically-plausible of the ETC to show that this simple model behaves the same as a fully mechanistic description of the ETC. Finally we conclude with suggestions for future research directions using simplified physically-plausible modelling as a strategy in systems and synthetic biology.

2 Modelling Bioenergetics of Biochemical Systems using Bond Graphs

Bond graphs provide a convenient modular framework for modelling energy flow within and across different physical domains: electrical, mechanical, chemical and so on; and as such are useful for representing biomolecular systems. In brief, bonds represent pairs of variables: potential and flow, whose product is power. In the biomolecular domain, the product of chemical potential $\mu$ (with units $\mathrm{J}\text{\,}{\mathrm{mol}}^{-1}$) and molar flow $v$ (with units $\mathrm{mol}\text{\,}{\mathrm{s}}^{-1}$) is power with units $\mathrm{J}\text{\,}{\mathrm{s}}^{-1}$ (Oster et al., 1971, 1973; Gawthrop and Crampin, 2014). Bonds connect components which represent either storage or dissipation of energy. In biomolecular systems, chemical potential is stored as concentration of chemical species, denoted **Ce ** 222In this paper, **Ce **components are used to represent chemical species and **C **components to represent electrical capacitors. , whereas chemical reactions, denoted **Re **, in which chemical species are converted from one form to another are dissipative processes. The biochemical network stoichiometry is represented in the coupling of **Ce ** components via the reactions **Re ** using bonds which represent the flow of energy, connected using common potential **0 **(‘zero’) and common flow **1 **(‘one’) junctions.

To illustrate, Figure 1 is the bond graph representation of the chemical reaction:

[TABLE]

**Ce ** components correspond to constitutive relations which relate the chemical potential to the amount of chemical species stored: the constitutive relations of Ce:$\mathbf{A}$ and Ce:$\mathbf{B}$ are:

[TABLE]

where $x_{A}$ and $x_{B}$ are the concentrations of \chA and \chB, $K_{A}$ and $K_{B}$ are species thermodynamic constants (${\mathrm{mol}}^{-1}$) for \chA and \chB, specific to each chemical species, $R$ is the universal gas constant and $T$ the absolute temperature.

The constitutive relations for the reaction components **Re ** provide the relationship between forward and backward chemical affinities $A$ (stoichiometric combinations of the chemical potentials) which provide the driving force for the reaction, and the molar flow (the reaction rate) $f$. The stoichiometry of reaction Re:$\mathbf{r}$ with formation of 2 molecules of species $B$ for each molecule of $A$ is represented by the two parallel bonds on the right hand side of Re:$\mathbf{r}$. With mass-action kinetics, the constitutive relation of Re:$\mathbf{r}$ is

[TABLE]

where $\kappa$ is a reaction rate constant ($\mathrm{mol}\text{\,}{\mathrm{s}}^{-1}$), specific to each reaction, and the forwards and backwards affinities are given by

[TABLE]

Combining these expressions gives the familiar mass-action expression for the reaction flow $f$:

[TABLE]

where the forward reaction rate constant $k^{+}=\kappa K_{A}$ and the reverse reaction rate constant $k^{-}=\kappa K_{B}^{2}$.

The bond graph approach is naturally allied to stoichiometric concepts (Klipp et al., 2016; Palsson, 2006, 2011, 2015). In particular, the stoichiometric matrix $N$ can be automatically generated from the network represented in the system bond graph. $N$ can be used to give species flows $f_{x}$ in terms of reaction flows $f$ and, conversely, reaction affinity $A$ in terms of species potentials $\mu$:

[TABLE]

In the case of the system of Figure 1:

[TABLE]

2.1 Modified mass action kinetics

Simplified representation of biomolecular system requires a representation of the reaction network that approximates, but does not fully represent the complete set of biochemical reactions. Physically-plausible models of biomolecular systems will therefore typically contain reactions which are an approximation to a sequence of elementary reactions. Thus even if elementary reaction steps have the mass action kinetics of Equation (2.4), this would not necessarily be the case for the overall reactions used to represent the system (see for example Atkins et al., 2018, chapter 17).

In bond graph terms, one may represent the dissipative reaction component with any appropriate constitutive relation for the reaction flow $f$ in terms of the forward and reverse affinities: mass action, as given by (2.4) leading to (2.7) is one example. In particular, non-elementary reactions may be represented using rate equations where, unlike the mass-action formulation, the concentration exponents are not the stoichiometric coefficients. One particular case of this would be to divide all of the stoichiometric coefficients by an positive integer constant $\alpha$ in the rate equations. Thus, for example, if $\alpha=2$, the reaction rate (2.7) corresponding to the reaction (2.1) would become:

[TABLE]

This corresponds to the non-integer stoichiometry

[TABLE]

Note that in the context of modelling the Mitochondrial Electron Transport Chain, the exponent $1/2$, corresponding to $\alpha=2$, commonly appears in the flux expression for complex III, as given by Beard (2005, equation B72) and Beard and Qian (2010, equation 7.38) for example, and the exponent $1/4$, corresponding to $\alpha=4$, appears in the flux expression for complex IV given by Beard (2005, equation B73) and Beard and Qian (2010, equation 7.41). This can be achieved by replacing the mass-action formula (2.4) by the modified mass-action (MMA) formula:

[TABLE]

which contains the additional parameter $\alpha$, which is used below as an essential part of the model fitting process.

For thermodynamic consistency, it is important that Equation (2.12) represents a dissipative system; that is, any non zero flow dissipates energy (Willems, 1972; Polderman and Willems, 1997; Willems, 2007). With this in mind, it is now shown that the MMA equation can be rewritten in mass-action form but with the positive constant $\kappa$ replaced by the positive function of concentration $\kappa_{\alpha}$. As a simple example of this, it can be verified that Equation (2.10) can be rewritten as

[TABLE]

As $x_{A}$ and $x_{B}$ are positive, $\kappa_{\alpha}$ is also positive. Thus (2.13) corresponds to the mass-action equation (2.4) with the positive constant $\kappa$ replaced by the positive function of concentration $\kappa_{\alpha}(x_{A},x_{B})$. The general modified-mass action kinetics of Equation (2.12) can also be rewritten in mass-action form with the positive constant $\kappa$ replaced by the positive $\kappa_{\alpha}$:

[TABLE]

2.2 Redox reactions

Oxidative phosphorylation involves a series of electrochemical redox reactions. Nicholls and Ferguson write that ‘Whereas all redox reactions can quite properly be described in thermodynamic terms by their Gibbs energy changes, electrochemical parameters can be employed because the reactions involve the transfer of electrons.” (Nicholls and Ferguson, 2013, chapter 3.3). Rather than to deal directly with conversion between electrical and chemical potentials and associated variables, it is convenient to have a common system of units and convert the chemical energy covariables chemical potential and molar flow to equivalent electrical energy covariables voltage and current (Gawthrop, 2017). The relevant conversion factor is Faraday’s constant $F\approx$96,485\text{,}\mathrm{C},\mathrm{m}\mathrm{o}\mathrm{l}^{-1}$$ (Nicholls and Ferguson, 2013; Gawthrop et al., 2017; Gawthrop, 2017). In particular, we define:

[TABLE]

Using these Faraday-equivalent variables, the **Ce **constitutive relations (2.2) and (2.3) become:

[TABLE]

and the modified mass-action formula (2.12) becomes:

[TABLE]

As noted by Nicholls and Ferguson, an advantage of transforming the chemical potentials into equivalent electrical potentials in the treatment of redox reactions is: “the ability to dissect the overall electron transfer into two half-reactions involving the donation and acceptance of electrons, respectively.” (Nicholls and Ferguson, 2013, chapter 3.3). For example, the two half-reactions:

[TABLE]

electron \che1- donation in the first (oxidation of A), and \che2- acceptance in the second (reduction of B), correspond to the overall reaction:

[TABLE]

Figure 2 shows a bond graph representation of these two half-reactions (LABEL:eq:Redox), explicitly representing the transfer of electrons \che- using the linear electrical capacitor represented by C:$\mathbf{E}$ with voltage $V$; the two-electron stoichiometry of reaction \chr1 is represented by the two parallel bonds.

If reaction \chr1 is in equilibrium, then the voltage $V$ is exactly that required to stop reaction \chr1 from proceeding and thus $V=-E_{1}$ where $E_{1}$ is the redox potential of reaction \chr1. Conversely, if reaction \chr2 is in equilibrium, then the voltage $V$ is exactly that require to stop reaction \chr2 from proceeding and thus $V=-E_{2}$ where $E_{2}$ is the redox potential of reaction \chr2.

2.3 Hierarchical Modelling

Hierarchical modelling and modularity provide one approach to understanding the complex systems associated with cellular biochemistry (Hartwell et al., 1999; Lauffenburger, 2000; Csete and Doyle, 2002; Bruggeman et al., 2002, 2008; Szallasi et al., 2010). Bond graphs provide an effective foundation for modular construction of hierarchical models of biochemical systems (Gawthrop et al., 2015; Gawthrop and Crampin, 2016). Bond graphs model the interaction between modules, in particular retroactivity (Jayanthi and Del Vecchio, 2011; Del Vecchio, 2013; Del Vecchio and Murray, 2014), in a straightforward manner whilst retaining thermodynamic compliance.

Bond graph modules use the notion of chemostats (Polettini and Esposito, 2014; Gawthrop and Crampin, 2016) which have an number of interpretations:

one or more species is fixed to give a constant concentration; this implies that an appropriate external flow is applied to balance the internal flow of the species. 2. 2.

as a **Ce **component with a fixed state. 3. 3.

as a module port through which chemical, mechanical or electrical energy flows.

Thus if the bond graph of Figure 2 were to be used as a module, then Ce:$\mathbf{A}$, Ce:$\mathbf{B}$, Ce:$\mathbf{C}$ and Ce:$\mathbf{D}$ could be chemostats. If the module were to be examined in isolation, then the interpretations of items 1 and 2 would be used; if, on the other hand, the module were to be embedded in a larger system, then the interpretation of item 3 would be used.

When examining the properties of a complex system, such as a whole cell model, the replacement of some modules by physically-plausible equivalents with the same ports would not only reduce computational complexity but also allow attention to be focussed on detailed models of other modules.

2.4 Dynamical Simulation

Bond graphs, together with the component constitutive relationships, can be used to automatically derive the ordinary differential equations (ODE) describing the system dynamics (Karnopp et al., 2012). These ODEs can be in symbolic form or in the form of computer code for a particular simulation engine. In general, a set of ODEs does not guarantee thermodynamic consistency; but, because these ODEs are derived from a bond graph, they inherit the thermodynamic properties of the bond graph.

In some systems, the system states are not independent; in particular, biomolecular systems usually have conserved moieties. In such circumstances, the system bond graph can be used to automatically generate the minimal number of ODEs describing the independent states from which dependent states and reaction flows can be derived (Gawthrop and Crampin, 2014, § 3(c)).

Although in this paper we have focused on species described by the standard logarithmic constitutive relation (2.2) and reaction flows determined by the mass action formula (2.4), bond graph components can have a wide range of constitutive relations constrained only by thermodynamics. A simple example of this is the modified mass-action formula of equation (2.12); a more complex example is the Goldman-Hodgkin-Katz flux equation used in bond graph models of action potential (Gawthrop et al., 2017). Moreover, models with complex characteristics, such as transporters, can be built from simple bond graph components (Pan et al., 2019).

2.5 BondGraphTools – a Python Toolkit

Computational tools are necessary for model capture, parameterisation and simulation. As the name suggests, BondGraphTools is an application programming interface (API) for capturing, simplifying and simulating bond graph models, and is an important part of the bond graph approach (Cudmore et al., 2019).

BondGraphTools is written in Python and is built upon the Scientific Python (SciPy) libraries, all of which are open source and easily accessible. The core use-case of BondGraphTools is to turn bond graphs into a set of reduced equations which can be then passed into other SciPy libraries (parameter estimation routines, or ODE integrators, for example). As model reduction is performed symbolically, the simplification routines are free from numerical errors, which is important for systems involving parameters that are unknown.

3 A Simplified Physically-Plausible Model for the Mitochondrial Electron Transport Chain

Mitochondria make use of redox reactions to provide the power driving many living systems. The key process in the generation of ATP is chemiosmotic energy transduction, whereby a sequence of redox reactions pumps protons across the mitochondrial inner membrane to generate the proton-motive force (PMF), an electrochemical gradient which is then used to power the synthesis of ATP. Generation of the PMF is accomplished by the mitochondrial electron transport chain. Beard and colleagues have developed the most comprehensive thermodynamically consistent models of mitochondrial oxidative phosphorylation including the electron transport chain (Beard, 2005; Wu et al., 2007; Bazil et al., 2016). Recently Gawthrop (2017) provided a bond graph model of mitochondrial oxidative phosphorylation based on the redox reactions associated with complexes \chCI, \chCIII and \chCIV of the mitochondrial electron transport chain.

In contrast, here we develop a simple, but physically-plausible, model based on the overall chemical reaction of the Electron Transport Chain in which \ch2 NADH is combined with \chO2 and \ch2 H+ to give \ch2 NAD+ and \ch2 H2O; the 2 protons (\ch2 H+) are consumed from the mitochondrial matrix. The free energy of this overall reaction pumps 20 protons across the mitochondrial inner membrane. Denoting the protons in the mitochondrial matrix as \chH_x+ and those in the mitochondrial inner membrane space as \chH_i+ the overall reaction is thus:

[TABLE]

Reaction (3.1) can be rewritten as the weighted sum of three reactions:

Table 2: Optimal parameters. The table summarises the fitting results using the two sets of simulation data ($\text{sim}_{0}$&$\text{sim}_{1}$) and the experimental data (exp) from Bazil et al. (2016). Rows 1–3 correspond to free $\alpha$ and the rows 4–6 to fixed $\alpha=6$.

Finally, the current-voltage relationship derived above and given in (LABEL:eq:explicit) is plotted in Figure 5 for different sets of fitted parameters given in Table 3, along with simulation results from the full model of Bazil et al. (2016), and the corresponding experimental data, showing that the explicit formula for the steady state current/voltage behaviour of the ETC is well captured by the physically plausible model. Figures 5(a)–5(c) correspond to rows 4–6 of Table 3; Figure 5(d) corresponds to row 3 of Table 3.

4 Discussion

Simplified models of biochemical and biophysical processes have a central role to play in the development of large whole-cell and multi-scale Physiome models, and the use of such models in biomedical and synthetic biology applications. Here we have argued that such simplified models need to be physically plausible, in the sense that they are consistent with the laws of physics (for example, that they obey mass conservation, are consistent with thermodynamic principles, and so on) as well as providing a suitable fit to available data sets. We have demonstrated that energy-based modelling using bond graphs provides a useful framework for the development of such models. The advantages of thermodynamically-consistent modelling have been argued recently by us and a number of other authors (Beard and Qian, 2010; Gawthrop and Crampin, 2014; Klipp et al., 2016; Gawthrop and Crampin, 2017; Mason and Covert, 2019). Key amongst these advantages are that models comprised of physically-plausible components are themselves physically plausible; that such models can be constructed and assembled in a modular fashion; and that such models enforce thermodynamic principles which allow identification of conserved moieties and furthermore restrict the possible parameter space. Here we have shown that by using a modified form of mass action we can generate a simple physically-plausible model of the mitochondrial electron transport chain that is able to reproduce experimentally measured properties of the system.

Bond graphs provide a framework within which is represented both the biochemical network stoichiometry and the constitutive relationship between thermodynamic driving force and biochemical reaction rate for each constituent bond graph element, describing the mechanism of enzymatic processes and so forth. Therefore, bond graphs describe the complete dynamical behaviour of the biochemical system, and can be used to derive the ordinary differential equations describing full system dynamics. A particular advantage, however, of a simple model as derived here over a fully mechanistic model is the possibility of deriving explicit formulae for key properties and behaviours of the system. Here we have shown that using the physically plausible modelling approach we can derive an explicit algebraic formula for the flux through the electron transport chain to the PMF that is generated across the mitochondrial membrane at steady state, given in equation (LABEL:eq:explicit). This is not in general possible from a full dynamical representation (whether represented as a bond graph or otherwise). This is of significance both as it drastically simplifies the model, and also because it allows a direct analysis of the dependence of the mitochondrial membrane potential on parameters of the ETC flux, as discussed below. Such a simplified but thermodynamically realistic representation of mitochondrial energy production is particularly suitable to be used in simulation of spatially-distributed networks of mitochondria (Jarosz et al., 2017; Ghosh et al., 2018), for which detailed mechanistic It remains to be determined how different possible modifications to mass action, or indeed other constitutive relations that may be used to relate chemical potential and reaction rate, affect the ability of simple physically plausible models to represent complex biochemical processes. models of mitochondrial bioenergetics have too high a computational overhead.

The physically plausible model of the electron transport chain derived above contains physical parameters which are known a-priori, as well as parameters which are model-dependent and therefore obtained by fitting to relevant data. In particular, the explicit formula (LABEL:eq:explicit) relating flow $f$ to potential difference $\Delta V$ has four known physical parameters: $\Delta V^{eq}$, $V_{N}$, $n_{1}$ and $n_{2}$ and four unknown parameters: $\alpha$, $K_{1}$, $K_{2}$ and $r_{loss}$. Using an optimization approach to model fitting, we have shown that the value of the parameter $K_{1}$ appears unimportant, as long as $K_{1}>10$; this corresponds to the rate constant $\kappa_{1}$ of reaction $r_{1}$ being large enough so that there is negligible potential drop across the reaction; in other words, this requires that the reaction (LABEL:eq:donate) is operating essentially at equilibrium under the experimental conditions considered.

In contrast, the $\alpha$ parameter of the modified mass-action equation (2.21) is found to be important. Choosing $\alpha=6$ gives a good fit for both the experimentally-fitted simulations and experimental data. This dependence is to be expected as the physically plausible model subsumes a number of individual reactions and this is known to lead to non-stoichiometric exponents (Atkins et al., 2018, chapter 17).

We have shown that despite it’s simplicity, the physically-plausible model fits the (complex) simulation data from Bazil et al. (2016) closely (RMS error $\epsilon\ll$1\text{,}\mathrm{mV}$$) for both values of \chPi and for free $\alpha$ and fixed $\alpha=6$. Furthermore the experimental data can be fitted with the model with an error $\epsilon$ of about $2\text{\,}\mathrm{mV}$.

Despite the success of the simplified, physically-plausible model of the electron transport chain that we have developed here, it should be noted that there are choices and trade-offs inherent in the simplification process. Therefore it is important that such simplified models are analysed and used only in the appropriate context. Firstly, the simplified model was generated using a specific form of modified mass action kinetics. While motivated by existing literature and approaches for representing non-elemental biochemical reaction steps, different choices could have been made. It remains to be determined how different possible modifications to mass action, or indeed other constitutive relations that may be used to relate chemical potential and reaction rate, affect the ability of simple physically plausible models to represent complex biochemical processes.

Secondly, the simplified model was developed using data relevant to normal physiological conditions. Most significantly, simplifying assumptions have been made about the physiological regime in which the model operates, and hence perturbations to species concentrations and specific enzymatic regulators outside of this regime are not captured in the simplified model. Full mechanistic exploration of mitochondrial dynamics under a broad range of perturbations and experimental data beyond those used to fit the simplified model should of course be pursued using the full bond graph representation of mitochondrial bioenergetics.

Finally, because the physically-plausible model is energy based, any or all of the **Ce **components can provide connections with energy-based models of other parts of the mitochondria system such as the TCA cycle, ATPase and ROS generation. In bond graph terminology, the **Ce **components become ports (Gawthrop and Crampin, 2018b) with which to connect to other bond graph components representing other aspects of mitochondrial biochemistry, for simulation and analysis of larger-scale models of mitochondrial and cellular bioenergetics. In our future work we intend to exploit this feature of bond graphs to investigate in larger models of mitochondrial function the basis of ROS generation and damage.

Acknowledgements

PJG would like to thank the Melbourne School of Engineering for its support via a Professorial Fellowship. This research was in part conducted and funded by the Australian Research Council Centre of Excellence in Convergent Bio-Nano Science and Technology (project number CE140100036). The authors would like to thank an anonymous reviewer for pointing out an error in an early version of the manuscript and suggesting a number of points of clarification.