A thermodynamic framework for modelling membrane transporters
Michael Pan, Peter J. Gawthrop, Kenneth Tran, Joseph Cursons, Edmund J. Crampin
AA thermodynamic framework for modelling membrane transporters
Michael Pan a , Peter J. Gawthrop a , Kenneth Tran b , Joseph Cursons c,d ,Edmund J. Crampin a,e,f, ∗ a Systems Biology Laboratory, School of Mathematics and Statistics, and Department of BiomedicalEngineering, Melbourne School of Engineering, University of Melbourne, Parkville, Victoria 3010, Australia b Auckland Bioengineering Institute, University of Auckland c Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, Parkville, Victoria 3052,Australia d Department of Medical Biology, School of Medicine, University of Melbourne, Parkville, Victoria 3010,Australia e School of Medicine, Faculty of Medicine, Dentistry and Health Sciences, University of Melbourne, Parkville,Victoria 3010 f ARC Centre of Excellence in Convergent Bio-Nano Science and Technology, Melbourne School ofEngineering, University of Melbourne, Parkville, Victoria 3010, Australia
Abstract
Membrane transporters contribute to the regulation of the internal environment of cells bytranslocating substrates across cell membranes. Like all physical systems, the behaviourof membrane transporters is constrained by the laws of thermodynamics. However, manymathematical models of transporters, especially those incorporated into whole-cell models,are not thermodynamically consistent, leading to unrealistic behaviour. In this paper we usea physics-based modelling framework, in which the transfer of energy is explicitly accountedfor, to develop thermodynamically consistent models of transporters. We then apply thismethodology to model two specific transporters: the cardiac sarcoplasmic/endoplasmic Ca ATPase (SERCA) and the cardiac Na + /K + ATPase.
Keywords: bond graph, biochemistry, chemical reaction network, biomedical engineering,systems biology
1. Introduction
The cell membrane is a physical barrier that separates the internal environment of acell from its external environment, as well as many compartments within the cell. Bytransporting chemical species across the membrane, the cell regulates concentration withineach compartment, providing the environment necessary for many cellular processes (Keener ∗ Corresponding author
Email addresses: [email protected] (Michael Pan), [email protected] (Peter J. Gawthrop), [email protected] (Kenneth Tran), [email protected] (Joseph Cursons), [email protected] (Edmund J. Crampin)
Preprint submitted to Elsevier Wednesday 13 th June, 2018 a r X i v : . [ q - b i o . B M ] J un nd Sneyd, 2009). A key contributor to these transport systems are membrane transporters,which are proteins located in cell membranes that maintain cell volume (Glitsch, 2001), setup ionic gradients required for electrical (Glitsch, 2001) and calcium signalling (Periasamyand Kalyanasundaram, 2007; Blaustein and Lederer, 1999), and transport sources of energyinto the cell (Mueckler and Thorens, 2013). As with all physical processes, transporters mustcomply with the principles of thermodynamics (Oster et al., 1973; Hwang, 2004; Beard andQian, 2008). Because each chemical species within a solution is associated with a chemicalpotential that increases with concentration, passive transporters can only move substratesfrom a region of high concentration to a region of low concentration (Keener and Sneyd,2009; Beard and Qian, 2008). To move a substrate against a concentration gradient, a sourceof energy must be provided. For example, active transporters use ATP hydrolysis to drivethe transport of substrates against a potential gradient (Keener and Sneyd, 2009).Mathematical models of membrane transporters have been developed for the purposeof understanding their transport mechanism, and predicting their behaviour beyond exper-imental measurements. However, despite the wealth of models available for transporters,thermodynamic consistency has not usually been applied. As a result, many of these modelsdescribe physically infeasible systems where for example, species are transported againsttheir potential gradients in the absence of an energy source, therefore generating energyout of nowhere. In the current literature, thermodynamic consistency is commonly violatedthrough the use of equations that describe irreversible transporters, neglect dependence oncertain metabolites, or have incorrect equilibrium points (Gawthrop and Crampin, 2014;Smith and Crampin, 2004; Tran et al., 2009). While there are methods for incorporatingthermodynamic constraints such as detailed balance (Liebermeister et al., 2010; Beard andQian, 2008; Smith and Crampin, 2004) and the Nernst potential (Keener and Sneyd, 2009),they tend to be scattered around the literature, and are not universally applied. Furthermore,many transporters are electrogenic, therefore there is an interaction between chemical andelectrical potential (Smith and Crampin, 2004; Terkildsen et al., 2007), and the multidomainnature of these transporters can confound efforts to develop thermodynamically consistentmodels. In some cases, thermodynamic inconsistency may impact on the ability of a modelto remain physiological under a wide range of conditions. For example, in the context ofheart failure, where ATP is depleted, active transporters such as the Na + /K + ATPase andsarcoplasmic/endoplasmic Ca ATPase (SERCA) operate at reduced rates. Whole cellmodels of cardiac cells generally do not describe this metabolite dependence because theyneglect thermodynamic constraints (Smith and Crampin, 2004), and while thermodynamicallyconsistent models capture these effects (Terkildsen et al., 2007; Tran et al., 2009), they tendto be the exception rather than the rule.To facilitate the incorporation of thermodynamic constraints into models of membranetransporters, we require a framework that is (a) thermodynamically consistent by design;(b) able to reproduce known thermodynamic constraints; and (c) general enough to model awide range of transporters. We propose the use of bond graphs, which are a physics-basedframework in which energy transfer is explicitly modelled through connections betweenphysical components, therefore satisfying thermodynamic consistency (Gawthrop and Smith,1996; Borutzky, 2010). The bond graph representation is also domain-independent, therefore2t is general enough to represent a wide range of physical systems. Bond graphs were originallyinvented by Henry Paynter for use in hydroelectric systems (Paynter, 1961), but they alsonaturally represent electrical and mechanical systems (Borutzky, 2010). The reader is referredto the texts by Gawthrop and Smith (1996), Borutzky (2010) and Gawthrop and Bevan(2007) for a comprehensive introduction to bond graph theory. More recently, bond graphshave been extended to chemical (Thoma and Bouamama, 2000), biochemical (Oster et al.,1973; Gawthrop and Crampin, 2014) and electrochemical systems (Gawthrop, 2017), enablingbond graph modelling of membrane transporters such as the sodium-glucose transport protein1 (SGLT1) (Gawthrop and Crampin, 2017). Models of electrogenic membrane transportershave recently been incorporated into a bond graph model of cardiac electrophysiology (Panet al., 2018a). A further advantage of modelling transporters as bond graphs is that they aremodular, therefore individual models of transporters are easily coupled to other bond graphmodels to build comprehensive models of biological systems (Gawthrop et al., 2015).In this paper, we review the bond graph theory required for modelling ion transportersand use bond graphs to build simple hypothetical models of transporters (Section 2), demon-strating that bond graphs capture important thermodynamic concepts such as detailedbalance, free energy of reaction and the Nernst potential. We then develop a bond graphmodel of cardiac SERCA based on existing work by Tran et al. (2009), and used the bondgraph model to assess energy consumption and efficiency (§ 3.1). Bond graphs also providea framework for detecting thermodynamic inconsistencies within existing models. In § 3.2we develop a new model of Na + /K + ATPase based on existing work by Terkildsen et al.(2007), and verify that the model complies with thermodynamic constraints. These examplesillustrate that bond graphs are a unifying framework for accounting for thermodynamicconstraints in models of membrane transporters. We believe that the bond graph approachwill prove to be a powerful tool in the development of thermodynamically consistent modelsof transporters, and other cellular processes in which energy transduction plays an importantrole.
2. Hypothetical models
In biochemical systems, Gibbs free energy (also called Gibbs energy or free energy) istransmitted between species and reactions (Atkins and De Paula, 2006). Like all physicalsystems, biochemical systems must comply with the laws of thermodynamics, thereforereactions can only proceed in the direction of decreasing chemical potential. However,chemical potential is ignored by many existing models, often resulting in physically unrealisticbehaviour (Gawthrop and Crampin, 2014; Gawthrop et al., 2015). Bond graphs accountfor chemical potential by decomposing the transfer of power into a product of chemicalpotential ( µ [J / mol]) and molar flow rate ( v [mol / s]) (Gawthrop and Crampin, 2014). Poweris transmitted between components within a bond graph by connecting them with bondsthat carry chemical potential and molar flow rate. Because the power leaving one componentis transmitted to another, bond graphs are thermodynamically consistent, that is, the flow ofenergy is explicitly accounted for. In bond graph terminology, chemical potential is known3 BC D E S i S e E E C:Si 0 0 00 11C:E1C:Se C:E2Re:r1Re:r2
Time (s) M o l a r a m oun t x E1 x E2 Time (s) -2024 R ea c t i on v e l o c i t y v v x Se -0.500.51 v cyc x Si = x Se Figure 1: A simple enzyme cycle. (A)
Kinetic model; (B)
Bond graph model. S i and S e are theintracellular and extracellular substrate respectively, E is the unbound transporter, and E is the transporterbound to substrate. A simulation was run with unity parameters, the initial conditions x E = x E = 1 and x Si = 100 and x Se = 10. The resulting molar amounts of transporter states (C) and reaction velocities (D) were plotted against time. (E) Using x Si = 100, the steady-state cycling rates were plotted against x Se , withthe blue dot indicating the equilibrium point.
4s the effort variable, and molar flow rate is known as the flow variable (Gawthrop andCrampin, 2014).In this section, we use a simple transporter model to outline the essential aspects ofmodelling biochemical systems using bond graphs (Figure 1A). This transporter modelcontains a two-state enzyme cycle, where an uncharged substrate S binds to the enzyme onthe intracellular (i) side, and unbinds on the extracellular (e) side. We denote the substratein a specific compartment as S x , where x ∈ { i , e } . The kinetic model is shown in Figure 1A,and the bond graph representation of this reaction cycle is shown in Figure 1B. Because athermodynamic approach implies reversibility for all reactions (Polettini and Esposito, 2014;Gawthrop and Crampin, 2014), the transport of the substrate can proceed in both directions,depending on the direction of the chemical gradient.Chemical energy is stored within the species of a biochemical system (Atkins and De Paula,2006). This energy allows the system to do work through chemical reactions. The ability of aspecies s to drive a reaction is given by the chemical potential µ s [J/mol], which is dependenton the molar amount of species x s [mol] in a logarithmic manner (Keener and Sneyd, 2009;Atkins and De Paula, 2006): µ s = RT ln( K s x s ) (1)where R = 8 .
314 J / K / mol is the ideal gas constant, T [K] is the absolute temperature and K s [mol − ] is the species thermodynamic constant. Therefore, the chemical potentials of thespecies in the enzyme cycle are µ Si = RT ln( K Si x Si ) (2) µ Se = RT ln( K Se x Se ) (3) µ E1 = RT ln( K E1 x E1 ) (4) µ E2 = RT ln( K E2 x E2 ) (5)Because species store energy, they are represented as C components in the bond graph model(Figure 1B), as they are analogous to capacitors in electrical circuits.The chemical energy stored within the various species of a biochemical system is usedto drive reactions, which convert chemical species into other chemical species, dissipatingchemical energy in the process. The rate of reaction can be related to the chemical potentialof its reactants and products using the Marcelin-de Donder equation (Gawthrop and Bevan,2007; Oster et al., 1973): v R = κ R ( e A fR /RT − e A rR /RT ) (6)where for a reaction R , v R [mol/s] is the rate of reaction, κ R [mol/s] is the reaction rateconstant, A fR [J/mol] is the forward affinity (representing the total chemical potential of thereactants) and A rR [J/mol] is the reverse affinity (representing the total chemical potential ofthe products). For the enzyme cycle example in Figure 1B, the reaction rates of the upper5nd lower reactions are given by the equations v = κ ( e A f /RT − e A r /RT ) (7) v = κ ( e A f /RT − e A r /RT ) (8)As seen in Figure 1B, reactions are represented as Re components in bond graph modelling.In the enzyme cycle, some species are involved in more than one reaction. Therefore werequire constraints to ensure that (a) the same chemical potential is used for the involvementof the same species in different reactions and (b) the contributions of each reaction aresummed when calculating the rate of change of each species. These constraints are capturedby the junction in bond graph modelling, which accounts for the first constraint by settingthe chemical potentials of all connected bonds to be equal. In Figure 1B, there are four junctions corresponding to each species. The second constraint is captured through theconservation equations v Si = − v (9) v Se = v (10) v E1 = v − v (11) v E2 = v − v (12) junctions are analogous to parallel connections in electric circuits, and the Eqs. 9–12represent the biochemical version of Kirchhoff’s current law. Note that the molar flow ratesof the bonds connecting to each junction sum to zero, therefore the junction conservespower (Borutzky, 2010).In this example, both reactions involve the combination or dissociation of chemical species,and therefore have multiple reactants or products. To account for this, we require constraintsso that (a) the rate at which each species is produced/consumed by the reaction is equalto the rate of reaction and (b) the affinities of each reaction are the totals of the chemicalpotentials of the reactants or products. These constraints are captured by junctions inbond graph modelling, which satisfy the first constraint by fixing the molar flow rates ofeach connected bond to equal values. In Figure 1B, there are two junctions representingthe association and dissociation of substrate. The second constraint is accounted for throughthe conservation laws A f = µ Si + µ E1 (13) A r = µ Se + µ E1 (14) junctions are analogous to series connections in electric circuits, and Eqs. 13–14 arebiochemical versions of Kirchhoff’s voltage law. Note that since the chemical potentials of allbonds connected to junctions sum to zero, they are power conserving, like the junction(Borutzky, 2010).Using the components described above, it is possible to derive a differential equation for6he enzyme cycle. By substituting Eqs. 13–14 into Eqs. 7–8, we find that the reaction ratesare v = κ ( e ( µ Si + µ E1 ) /RT − e µ E2 /RT ) = κ K Si K E1 x Si x E1 − κ K E2 x E2 (15) v = κ ( e µ E2 /RT − e ( µ Se + µ E1 ) /RT ) = κ K E2 x E2 − κ K Se K E1 x Se x E1 (16)Hence the bond graph components are able to reproduce the mass action equations, with thekinetic parameters k +1 = κ K E K Si (17) k − = κ K E (18) k +2 = κ K E (19) k − = κ K E K Se (20)By using the conservation laws from Eqs. 9–12, the rates of change for each species are˙ x Si = v Si = − v = − k +1 x Si x E1 + k − x E2 (21)˙ x Se = v Se = v = k +2 x E2 − k − x Se x E1 (22)˙ x E1 = v E1 = v − v = − k +1 x Si x E1 + k − x E2 + k +2 x E2 − k − x Se x E1 (23)˙ x E2 = v E2 = v − v = k +1 x Si x E1 − k − x E2 − k +2 x E2 + k − x Se x E1 (24)The above equations are applicable to an isolated system. However, since transportersare dissipative systems, i.e. they release energy as heat, it is impossible for these to operatecontinuously at steady state without an external power supply. Therefore, we model theseexternal power supplies by holding certain species at constant concentrations rather thanallowing them to dynamically vary. Such species are known as chemostats (Polettini andEsposito, 2014). In this example, we assume that the concentrations of the substrates S i andS e are constant, so that Eqs. 21 and 22 are replaced by˙ x Si = 0 (25)˙ x Se = 0 (26)Because an external flow of species is required to keep the concentration of the speciesconstant, chemostats are not energy conserving, but rather represent the influx or release ofenergy into the external environment. Thus chemostats turn a biochemical system from anisolated system into an open system that communicates with its external environment, andcan therefore be used to couple models together. In bond graph modelling, C components7re replaced with Se (effort source) components when they become chemostats to indicatethe transfer of power with an external source. Since species can be represented using C or Se components interchangeably depending on the purpose of the model (Gawthrop et al., 2017;Gawthrop, 2017), we represent species as C components in all diagrams within this paper,and mention in text whether certain species are treated as “chemostats” in the analysis.For the remainder of Section 2, we assume that:1. Apart from the states of the transporter (e.g. E , E ), the concentrations of all speciesare constant (i.e. they are modelled as chemostats).2. The volumes of all compartments are equal, so that the amount of each species directlycorresponds to concentration. In real biological systems, compartments will havedifferent volumes, and we show how to incorporate these effects in Appendix B of theSupplementary Material.3. The parameters K and κ for each species and reaction take on unity values.Assumptions 2 and 3 are made to simplify the analysis in this section, although we showhow to generalise beyond these assumptions in later sections.Figure 1C,D shows simulations of this model when the amount of enzyme is small relativeto the amount of substrate, thus we expect the system to achieve a steady state relativelyquickly. At steady state, the amounts of each pump state are constant, and the two reactionvelocities converge towards the same value, as dictated by the pathway analysis of Gawthropand Crampin (2017). We are often interested in the cycling (or turnover) rate at steadystate, v cyc , which is given by v cyc = V /e , where V is the steady-state reaction rate, and e = x E1 + x E2 is the total amount of transporter (Atkins and De Paula, 2006).Figure 1E shows the effect of changing x Se on cycling rate, indicating an inverse relation-ship: as x Se increases, the transporter operates at a lower rate, and eventually operates inthe reverse direction as expected. The bond graph model captures the fundamental physicalconstraint that the equilibrium point between the forward and reverse regimes of operationoccurs when x Se = x Si . Therefore the simple transporter can only allow transport of itssubstrate down its concentration gradient.The direction in which each reaction proceeds is determined by the Gibbs free energyof reaction ∆ G . Because reactions can only run in the direction of decreasing chemicalpotential, they will only proceed in the forward direction if ∆ G is negative. The Gibbs freeenergy of reaction relates to the affinities A = A f − A r of the Re components in bond graphmodels. Since efforts in the biochemical domain are associated with Gibbs free energy, theGibbs free energy of reaction is the negative of affinity:∆ G = − A = A r − A f (27)Therefore, the free energies of each of the reactions in this example are∆ G = µ E2 − µ E1 − µ Si (28)∆ G = µ E1 + µ Se − µ E2 (29)8t steady state, the transporter’s direction of operation is determined by the Gibbs freeenergy of the overall reaction S i (cid:10) S e , which is also the sum of all reactions in its cycle(Gawthrop and Crampin, 2017):∆ G = ∆ G + ∆ G = µ Se − µ Si (30)By substituting Eq. 1 and setting ∆ G = 0, the equilibrium of the system can be found:∆ G = RT ln( K Se x Se ) − RT ln( K Si x Si ) = RT ln ( x Se /x Si ) = 0 (31) ⇒ x Se = x Si (32)Therefore, as expected, the point at which the free energy is zero corresponds to theequilibrium of the transporter.To ensure that the equilibrium occurs at x Se = x Si , we must specify the equilibriumbetween the substrate in each compartment: K Si /K Se = K eq = 1 (33)which is exactly equivalent to the detailed balance constraint used in kinetic models (Keenerand Sneyd, 2009; Smith and Crampin, 2004; Tran et al., 2009): k +1 k +2 k − k − = 1 (34)(this can be easily checked by substituting Eq. 33 into Eqs. 17–20). Because bond graphs arethermodynamically consistent, the thermodynamic constraint is simpler and more intuitivewhen compared to that for the kinetic parameters. In the previous example, we observed that a passive transporter was unable to move asubstrate against its concentration gradient. However, many transporters, including activetransporters and exchangers, are able to move substrates against their concentration gradients.Such transporters couple the movement of substrate against a concentration gradient to aprocess that provides sufficient energy to enable transport. In Figure 2A, we show a simplemechanism for coupling together these processes. In addition to binding S i and unbinding S e ,the transport also binds A and unbinds B , giving rise to the overall reactionS i + A (cid:10) S e + B (35)The reaction A (cid:10) B provides energy for the transport of substrate, and may for examplerepresent ATP hydrolysis in active transporters, or the transport of another species down itsconcentration gradient in exchangers or cotransporters. The bond graph representation ofthis chemical reaction network is given in Figure 2B.To achieve the correct equilibrium point, we use the following thermodynamic constraint,9 E S i E E AS e B ABC D E -1 -0.5 0 0.5 1 G -0.01-0.00500.0050.01 v cyc G = 0, v ss = 0 Varying AVarying BVarying Se x A -0.04-0.020 v cyc x A = 10 x Se -0.0500.050.1 v cyc x A = 2x A = 5x A = 10x A = 20 Figure 2: Transport of a species by coupling to an energy supply. (A)
Kinetic model; (B)
Bondgraph model. S i and S e are the intracellular and extracellular substrate respectively, A and B are speciesthat provide power for the transport cycle, and E , E , E and E are the transporter states. Simulationswere run using the initial conditions x E = x E = x E = x E = 0 .
5, and chemostats x B = 1, x Si = 10. (C) The effect of x Se on steady state cycling rate for different x A . Equilibrium points, where the cycling ratesare zero, are denoted by the dots. (D) Effect of x A on the cycling rate, with x Se = 100. The equilibrium(blue dots) occurs at x A = x Se / ( x Si x B ) = 10. (E) The relationship between the Gibbs free energy of thetransporter and cycling rate, when A, B and S e are varied. K Si /K Se = 1 (36) K A /K B = 1 (37)Note that K A /K B can be set to different values depending on the process that A (cid:10) B represents; it has a value of 1 if it represents the transport of an uncharged substrate (asthe equilibrium point occurs when the substrate has the same concentration on either sideof the membrane), whereas it is determined by the standard Gibbs free energy of reactionif it represents a reaction such as ATP hydrolysis (see Appendix B of the SupplementaryMaterial).For a kinetic model with the parameters k +1 , k − , k +2 , k − , k +3 , k − , k +4 and k − , thecorresponding thermodynamic constraint is k +1 k +2 k +3 k +4 k − k − k − k − = K Si K Se K A K B = 1 (38)We note here that the thermodynamic constraint between kinetic parameters is more compli-cated than that in § 2.1 because there are more reactions in the cycle. When the numberof reactions in biochemical cycles increases, these thermodynamic constraints become morecomplex, and harder to derive. By contrast, when bond graph parameters are used, thermo-dynamic constraints remain simple regardless of cycle complexity, and are more intuitive toformulate.In Figure 2C we simulated the bond graph model to show the steady-state cycling ratesunder different values of x A and x Se . We note that for a passive transporter the equilibriumpoint occurs at x Se = x Si = 10, however in this example with x A > x B , the equilibrium points(represented by dots) occur at x Se > x Si . With all thermodynamic parameters set to 1, theequilibrium is given by the equation x Se = x Si x A /x B (39)Thus when x A > x B (i.e. the coupled process releases energy), x Se is greater than x Si atequilibrium. Therefore there is a region of operation (when x Si < x Se < x Si x A /x B ) where thetransporter can transport S against its concentration gradient. As the value of x A increases,the equilibrium point shifts towards the right (Figure 2C), indicating that the transporteris able to move S against a greater concentration gradient. This is because with a greater x A , the reaction A (cid:10) B provides more power, driving the transport in the forward direction(Figure 2D). 11he Gibbs free energies of each of the reactions are∆ G = µ E2 − µ E1 − µ Si (40)∆ G = µ E3 − µ E2 − µ A (41)∆ G = µ E4 + µ Se − µ E3 (42)∆ G = µ E1 + µ B − µ E4 (43)and the free energy of the transporter is∆ G = ∆ G + ∆ G + ∆ G + ∆ G = µ Se + µ B − µ Si − µ A (44)By substituting Eq. 1 and setting ∆ G = 0, we can recover the equilibrium relationship inEq. 39: ∆ G = RT ln( K Se x Se ) + RT ln( K B x B ) − RT ln( K Si x Si ) − RT ln( K A x A )= RT ln (cid:18) K Se x Se K B x B K Si x Si K A x A (cid:19) = RT ln (cid:18) x Se x B x Si x A (cid:19) = 0 (45) ⇒ x Se x B x Si x A = 1 (46)Therefore, the equilibrium for steady-state operation corresponds to the point where theGibbs free energy of the transport process is zero. We verify this fundamental physicalconstraint by plotting Gibbs free energy against cycling rate (Figure 2E). We varied theconcentrations of A, B and S e to generate three different curves. Despite differences in theshape of each curve, they each pass through the equilibrium point ∆ G = 0, v cyc = 0, verifyingthat the bond graph model correctly captures the relationship between Gibbs free energyand equilibrium. Furthermore, the transporter operates in the forward direction ( v cyc > v cyc <
0) with apositive (unfavourable) free energy.
Many membrane transporters, including ion exchangers, cotransporters and active iontransporters, move a charged species across a membrane. For those membranes within thecell (such as the plasma membrane, and inner mitochondrial membrane) that are charged,the membrane potential influences both the kinetics and thermodynamics of transport. Inthis section we explore the impact on transport by modifying the enzyme example in § 2.1,so that the transported substrate has a single unit of positive charge (i.e. z = 1) (Figure 3A).For simplicity, we choose to assign the entirety of the electrical dependence to the forwardside of the first reaction, which results in an exponential dependence on membrane voltage,arising from thermodynamic constraints (Keener and Sneyd, 2009; Smith and Crampin,2004). 12 V m -1-0.500.51 v cyc A B S i S e E E C D E
Re:r100 00 01C:mem C:E2Re:r2C:Se C:E1C:Si 1TF:zF -0.1 -0.05 0 0.05 0.1 V m -500005000100001500020000 G -1 0 1 2 G -1-0.500.51 v cyc G = 0, v ss = 0 Figure 3: Transport of a charged species (A)
Kinetic model; (B)
Bond graph model. S +i and S +e are the intracellular and extracellular substrate respectively, E is the unbound transporter, and E isthe transporter bound to substrate. Simulations were run with the initial conditions x E = x E = 1, andchemostats x Si = 10, x Se = 100. (C) Relationship between the membrane potential and cycling rate. Theblue dot indicates the equilibrium potential, as predicted by the Nernst equation. (D)
The effect of membranepotential on the Gibbs free energy of the transporter. (E)
Relationship between the Gibbs free energy andcycling rate under the simulation conditions, with varying membrane potential. C component representing membrane capacitance.This C component has similar energy-storage properties to the C component for each chemicalspecies, however, in place of Eq. 1, it has a linear constitutive equation V m = q m /C , where V m [V] is the membrane potential and q m [C] is the charge. We use a membrane capacitanceof C = F . Since the capacitor is an electrical component, its power must be converted tochemical power to appropriately describe its influence on reaction kinetics. This conversionis given by Faraday’s constant F = 96485 C / mol, which relates the membrane potential andcurrent to chemical potential and molar flow rate (Gawthrop et al., 2017): µ m = zF V m (47) I = zF v (48)These transformations are described by the TF component in Figure 3B. Since µ m v = V m I ,the TF component is a power-conserving transformation.By substituting the relevant constitutive equations, the rate of reaction 1 is given by v = κ ( e ( µ Si + µ E + µ m ) /RT − e µ E /RT ) = κ K Si K E x Si x E e F V m /RT − κ K E x E (49)which corresponds to the exponential dependence of the kinetic reaction scheme. Becauseadding an electrical component affects the thermodynamics of transport, the equilibriumbecomes dependent on membrane potential. As a result, the kinetics of the transporter mustdepend on membrane voltage to account for changes in equilibrium, although in practicemodellers often make assumptions about where the electrical dependence lies due to a lackof experimental data (Smith and Crampin, 2004).Bond graphs incorporate thermodynamic constraints for electrogenic transport. Becausethe chemical reaction structure is similar to the system in § 2.1, the constraint in Eq. 33holds. However, because a single unit of charge is moved across the membrane in a singlecycle, there is an additional constraint for charge: z f − z r + z f − z r = z = 1 (50)where z indicates charge, subscripts indicate the reaction number, and superscripts indicatethe side of the reaction. In this example, z f = 1 and z r = z f = z r = 0, therefore Eq. 50 issatisfied.At equilibrium, the free energy of the overall reaction is zero (Gawthrop and Crampin,2017), and because K Si = K Se (Eq. 33),∆ G = µ Se − µ Si − µ m = RT ln( K Se x Se ) − RT ln( K Si x Si ) − F V m = 0 (51)14 Ca E i H ATPE i H ATPE i HCa ATP-E i Ca HP (cid:1) E sr Ca HP (cid:0) E sr HP (cid:2) E sr H HP (cid:3) E sr MgADPP P P P P P P P P R1-2 R2-4 R4-5R6-8R8-9R9-10 HATPE i HMgATP nH + nH + H + Pi H + R5-6R10-1 A C:MgATPC:P1 C:P2 C:P4 C:P6C:P8 C:P5C:P10 R e : R - R e : R - R e : R - R e : R - R e : R - R e : R - C:P9Re:R10-10 1 10 0 1 001010101 C:CasrRe:R2-2a10C:P2a 0 C:HTF:nCaiC:Cai TF : n Re:R5-6 C:MgADP1TF:nCasrC:Pi TF:n B Figure 4: Kinetic and bond graph representations of the cardiac SERCA model from Tranet al. (2009) . (A) Kinetic model, adapted from Fig. 2 of Tran et al. (2009). The numbers for each pumpstate are shown in blue boxes, and the names for each reaction are labelled in green. (B)
Bond graphrepresentation, with pump states shown in blue, and reactions in green. the membrane potential at equilibrium is V m = RTF ln (cid:18) x Se x Si (cid:19) (52)which is the familiar Nernst potential (Keener and Sneyd, 2009).In Figure 3C, this system was simulated under varying membrane potentials. Becausethe species is positively charged, changing the membrane potential in the positive directiondrives transport of the species out of the cell. The transporter achieved an equilibrium at amembrane potential of approximately 0.062, consistent with the Nernst equation in Eq. 52.In Figure 3D we plot the Gibbs free energy (defined in Eq. 51) against membrane potential.The membrane potential has a linear contribution to the free energy of the transporter,becoming more negative for positive membrane potentials due to the positively chargedsubstrate. Importantly, the zero of the Gibbs free energy coincides with the zero of thecycling rate in Figure 3C, thus the transporter is dissipative, and therefore thermodynamicallyconsistent. Figure 3E shows this directly, where the transporter operates in the positivedirection only if the Gibbs free energy is negative.
3. Thermodynamic models of cardiac cell transporters
The sarcoplasmic/endoplasmic Ca ATPase (SERCA) is an active ion transporter thatpumps calcium from the cytosol into the sarcoplasmic reticulum (SR), an intracellular Ca ] sr > [Ca ] i ), it couples calciumtransport to the hydrolysis of ATP in order to obtain the energy required for transport. Theoverall reaction of the pump is2Ca + MgATP (cid:10) + MgADP + P i + H + (53)Since the SR membrane is uncharged, the pump is driven purely by chemical energy.SERCA has been included in a wide variety of models of cardiac cell Ca cycling,excitation-contraction coupling and electrophysiology. Our analysis here is based on themodel by Tran et al. (2009) which has since been incorporated into a number of subsequentcardiac cell models because it describes the dependence of cycling rate on all metabolites(Tran et al., 2015; Williams et al., 2011; Walker et al., 2014). The chemical reaction networkand bond graph representation of the Tran et al. (2009) model are shown in Figure 4. Thestructure of the bond graph closely resembles the cyclic nature of the chemical reactionscheme. However, the bond graph representation explicitly shows that the H + ions involved inmultiple reactions are linked. Note that the bond graph representation uses TF componentsto describe multiple copies of a species in a single reaction; this is discussed further inAppendix A of the Supplementary Material.Discussion of how bond graph parameters may be determined from an existing kineticmodel are presented in Appendix B of the Supplementary Material. The kinetic parameterswere thermodynamically constrained, therefore it was possible to find a set of correspondingbond graph parameters by using Eq. B.9 of the Supplementary Material, with the stoichio-metric matrix and resulting bond graph parameters given in Appendix D. A comparison ofthe kinetic and bond graph models (Figure 5A) confirms that the two models match closely,with only minor discrepancies at higher cycling rates due the assumption of rapid equilibriumin the Tran et al. (2009) model.SERCA accounts for approximately 10–15% of energy consumption within cardiomyocytes(Tran et al., 2015; Schramm et al., 1994), and is seen as a major contributor to myocardialenergy expenditure. In heart failure, SERCA activity decreases, which causes a higherproportion of Ca to be removed via an alternative pathway – the less efficient Na + /Ca exchanger (NCX) – increasing energy expenditure (Kawase and Hajjar, 2008). Thereforeit is important to introduce the notion of energy and efficiency into models of transportersin cardiac cells to incorporate the metabolite dependencies required for studying changesin flux under heart failure, and to compare the efficiencies of transporters that move thesame substrates. While the parameters for kinetic models can be chosen such that they arethermodynamically consistent, energy-related quantities such as Gibbs free energy and powerconsumption do not arise naturally from these parameters. Because bond graphs explicitlymodel energy transfer, an advantage of using this framework is that the power consumptionand efficiency are easily calculated from the model. At steady state, the Gibbs free energyof the pump can be calculated from the chemical potentials of metabolites (Gawthrop and16 BC D [Ca ] sr (mM) -20-1001020 G ( kJ / m o l ) [Ca ] sr (mM) P o w e r( kJ / m o l / s ) [Ca ] sr (mM) E ff i c i en cy ( % ) [Ca ] sr (mM) C yc li ng r a t e ( s - ) KineticBond graph
Figure 5: Simulation of the SERCA pump. (A)
Comparison of cycling rates for kinetic and bondgraph models, reproducing part of Fig. 13 in Tran et al. (2009); (B)
Gibbs free energy; (C)
Powerconsumption per mol of pump; (D)
Pump efficiency. Simulations were run with [Ca ] i = 150 nM, pH = 4,[MgADP] = 0 . . G = 2 µ Casr + µ MgADP + µ Pi + µ H − µ Cai − µ MgATP (54)The relationship between free energy and SR Ca is shown in Figure 5B. Note that theGibbs free energy is zero at the equilibrium point of the pump in Figure 5A, as expectedof a thermodynamically consistent system. It is important to note that when ∆ G > flows from the SR to cytosol, and ATP is synthesised (Makinose and Hasselbach, 1971).We observe in Figure 5B that a thermodynamic framework captures this reversal mode ofoperation, as it is a fundamental physical constraint. The product of free energy and reactionrate gives the rate of power dissipation. As seen in Figure 5C, the power consumption ispositive under all conditions except at equilibrium, where it is zero.Because the bond graph approach can split energetic contributions from different sources,it is possible to assess the efficiency of the pump. For this, we define the affinity of ATPhydrolysis as A hyd = µ MgATP − µ MgADP − µ Pi − µ H (55)and the affinity of calcium transport as A tr = 2 µ Cai − µ Casr (56)Using the notion of pumping efficiency introduced in Gawthrop and Crampin (2018), wedefine the efficiency ρ as ρ = − A tr /A hyd , A hyd ≥ − A tr − A hyd /A tr , A hyd < − A tr (57)Thus if the pump operates in the forward direction, efficiency is the proportion of energyfrom ATP hydrolysis that is converted into energy for calcium transport, and if the pumpoperates in the reverse direction, efficiency is the proportion of energy supplied from calciumtransport that is used to generate ATP. We therefore expect a cusp point at A hyd = − A tr .The efficiency of the pump is plotted in Figure 5D. This SERCA model operates reasonablyefficiently under the simulated conditions, ranging from 70–100%, consistent with previousestimates that the pump is 85-90% efficient (Pinz et al., 2011). We note also a negativerelationship between cycling rate and efficiency, with the pump becoming less efficient thefurther it moves away from equilibrium. Conversely, the pump approaches an efficiency of100% as it nears equilibrium. However, it should be noted that in reality, SERCA pumps mayexhibit “slippage”, where the pump hydrolyses ATP without transporting any Ca (de Meis,2002). Incorporating such behaviour into the model would likely reduce the maximumoperating efficiency of the pump (Gawthrop and Crampin, 2017).18 TP (cid:1) E i K ATP (cid:0) E i K ATP (cid:2) E i ATP (cid:3) E i Na ATP (cid:4) E i Na ATP (cid:5) E i Na P (cid:6) E i (Na )P (cid:7) E e Na P (cid:8) E e Na P (cid:9) E e Na P (cid:10) E e K P (cid:11) E e K E e (K )ATP (cid:12) E e (K ) P (cid:13) E e K + K + K + K + Na + Na + Na + Na + Na + Na + MgATP P P P P P P P P P P P P P P P R1 R2 R3 R4 R5 R6R7R8R9R10R11R12R13R14R15
C:Ki C:NaiC:P1 C:P2 C:P3 C:MgADPC:P7C:P9C:P10C:P11C:H C:P_iC:P14C:MgATPC:P15 C:P4 C:P5 C:P6C:P8C:P13 R e : R R e : R R e : R R e : R R e : R R e : R R e : R R e : R R e : R R e : R C:P12 Re:R6Re:R7Re:R13Re:R14Re:R150 10 10 0 1 0 1 0 1 1 0100110 0101010101010 TF:zF_5TF:zF_80 0 00 0C:memC:Ke C:Nae AB Figure 6: The cardiac Na + /K + ATPase model. (A)
Kinetic model, with numbers for each pumpstate ( blue boxes ) and reaction names ( green ) are labelled, with corrected parameters shown in red. (B)
Bond graph model, with pump states coloured in blue, and reactions coloured in green.
150 -100 -50 0 50 100
Membrane voltage (mV) C u rr en t D en s i t y ( A / F ) [Na e ] = 50mM[Na e ] = 100mM[Na e ] = 150mM -150 -100 -50 0 50 100 Membrane voltage (mV) C yc li ng v e l o c i t y ( s - ) [Na e ] = 1.5mM[Na e ] = 50mM[Na e ] = 100mM[Na e ] = 150mM A B
Figure 7: Fit of the cardiac Na + /K + ATPase model to current-voltage measurements. (A)
Comparison of the model to extracellular sodium and voltage data (Nakao and Gadsby, 1989, Fig. 3), withcycling velocities scaled to a value of 55s − at V = 40mV. (B) Comparison of the model to whole-cell currentmeasurements (Nakao and Gadsby, 1989, Fig. 2A). [Na + ] i = 50 mM, [K + ] i = 0 mM, [K + ] e = 5 . .
4, [Pi] tot = 0 mM, [MgATP] = 10 mM, [MgADP] = 0 mM, T = 310 K. + /K + ATPase
The Na + /K + ATPase is responsible for maintaining the Na + and K + gradients that driveionic currents during the action potential of cardiac cells and many other excitable cells.Na + and K + are pumped against their electrochemical gradients, therefore their transportrequires the supply of energy from ATP hydrolysis. The overall reaction of the pump is3Na +i + 2K +e + MgATP (cid:10) +e + 2K +i + MgADP + P i + H + (58)In contrast to SERCA, in which only chemical potentials determine the direction of operation,the Na + /K + ATPase is driven by both chemical and electrical potentials because the plasmamembrane is electrically charged. Therefore a thermodynamically consistent model of Na + /K + ATPase must account for both the free energy of ATP hydrolysis as well as the energeticcontribution of the membrane potential.In this section, we outline a new model of cardiac Na + /K + ATPase, based on the earliermodel by Terkildsen et al. (2007) (Figure 6). While the model of Terkildsen et al. (2007) isa biophysically detailed model that incorporates some thermodynamic principles, the finalmodel is thermodynamically inconsistent. We updated the model, correcting the equilibriumconstants used for identical binding sites, the equilibrium constant for ATP hydrolysis andmathematical expressions arising from the rapid equilibrium approximation, as described inAppendix C of the Supplementary Material.We refitted the new model to the data used to parameterise the original Terkildsenet al. (2007) model (see Appendix C of the Supplementary Material for further detail).Comparisons of model simulations to data are given in Figures 7,8. The updated modelprovides a a good fit to the data, and has a comparable fit compared to the original model.20 B [MgATP] (mM) R e l a t i v e cyc li ng v e l o c i t y C [Na + ] i (mM) R e l a t i v e cyc li ng v e l o c i t y [K + ] e (mM) R e l a t i v e cyc li ng v e l o c i t y Figure 8: Fit of the cardiac Na + /K + ATPase model to metabolite dependence data.
Simulationconditions are displayed on the right of each figure. (A)
Comparison of the model to data with varyingintracellular sodium concentrations (Hansen et al., 2002, Fig. 7A), normalised to the cycling velocity at[Na + ] i = 50 mM. (B) Comparison of the model to data with varying extracellular potassium (Nakao andGadsby, 1989, Fig. 11A), normalised to the cycling velocity at [K + ] e = 5 . (C) Comparison of themodel to data with varying ATP (Friedrich et al., 1996, Fig. 3B), normalised to the cycling velocity at[MgATP] = 0 .
10 0 10 20
G (kJ/mol) -0.4-0.200.20.40.6 C yc li ng r a t e ( s - ) A B -400 -350 -300 -250 -200 V m (mV) -0.500.51 C yc li ng r a t e ( s - ) [MgATP]=1mM[MgATP]=6.95mM Figure 9: Simulation of the Na + /K + ATPase. (A)
Cycling rates of the pump near reversal potential; (B)
Relationship between Gibbs free energy and cycling rate. The curves represent different concentrationsof MgATP, from a concentration of 1 mM on the right, with increments of 1 mM up to a concentrationof 5mM on the left. The Gibbs free energy was varied by changing the membrane potential. For (A) and(B), simulations were run using [Na + ] i = 10 mM, [Na + ] e = 140 mM, [K + ] i = 145 mM, [K + ] e = 5 . . . .
95 mM, [MgADP] = 0 .
035 mM. Each pump state wasinitialised to 1/15 fmol, and steady states were estimated by running each simulation to steady state.
The fit in Figure 7A was slightly worse than the original model at lower extracellular Na + concentrations, but the model seems to be more consistent with experimental evidence thatthe saturated pump velocity has little sensitivity to extracellular Na + at positive membranepotentials (Nakao and Gadsby, 1989).The updated model is thermodynamically consistent, therefore it has a bond graphrepresentation (Figure 6B). The definition of bond graph parameters and stoichiometricmatrices from the fitting process are given in Appendix E of the Supplementary Material.At steady state, the Gibbs free energy of the cycle is∆ G = 3 µ Nae + 2 µ Ki + µ MgADP + µ Pi + µ H − µ Nai − µ Ke − µ MgATP − µ mem (59)= ∆ G + RT ln [MgADP][Pi][H + ][MgATP] ! + 3 RT ln [Na + ] e [Na + ] i ! + 2 RT ln [K + ] i [K + ] e ! − F V m (60)To verify the thermodynamic consistency of the bond graph model, we simulated the steady-state operation of the pump at voltages near equilibrium, under physiological (6.95 mM) andischaemic (1 mM) MgATP concentrations (Figure 9A). The equilibrium potentials of thepump, where cycling rates are zero, are indicated by dots. An expression for the equilibrium22otential can be derived by solving for ∆ G = 0: V eq m = 1 F ∆ G + RTF ln [MgADP][Pi][H + ][MgATP] + 3 ln [Na + ] e [Na + ] i + 2 ln [K + ] i [K + ] e ! (61)The predicted equilibrium potentials of − . − . v cyc = 0 when ∆ G = 0, and v cyc is positive only at negative free energies.Thus Figure 9 verifies that the updated model is thermodynamically consistent. Like SERCA,the Na + /K + ATPase is reversible, and has been experimentally observed to synthesise ATPunder artificially large ionic gradients (Glynn, 2002). A thermodynamic framework is idealfor describing such phenomena.
4. Discussion
In this paper we have discussed the essential thermodynamic principles underlyingmembrane transporters, and how the bond graph framework captures thermodynamicconstraints in these systems. We illustrated using hypothetical models that bond graphscan be used to model simple membrane transport mechanisms, incorporating relevantphysical and thermodynamic constraints for both nonelectrogenic and electrogenic transport.The bond graph approach provides a single framework that naturally incorporates knownthermodynamic constraints in these systems, including detailed balance (Keener and Sneyd,2009; Liebermeister et al., 2010) and the Nernst potential (Keener and Sneyd, 2009).We applied the bond graph approach to two membrane transporters. Using the model ofcardiac SERCA by Tran et al. (2009), we demonstrated that thermodynamically consistentmodels can be represented as bond graphs, and that energetic quantities are easily calculatedfrom the bond graph model (§ 3.1). We then developed a bond graph model of cardiacNa + /K + ATPase based on earlier work by Terkildsen et al. (2007) and Smith and Crampin(2004).
Because bond graphs are a general-purpose modelling framework for physical systems,they unify a number of known thermodynamic constraints in the literature under a singleframework, allowing thermodynamic constraints to be applied more consistently across modelsof biological systems. Whereas in kinetic models, Wegscheider conditions and detailed balanceare used to constrain kinetic parameters in models of membrane transporters (Keener andSneyd, 2009; Smith and Crampin, 2004; Tran et al., 2009), the thermodynamic parameters ina bond graph model automatically account for this constraint (Gawthrop and Crampin, 2014).In electrogenic transport, the Nernst equation (Keener and Sneyd, 2009) is a thermodynamicconstraint that captures the electrochemical equilibrium. This constraint is automatically23ccounted for in bond graph modelling, because the Nernst potential is automatically derivedfrom the constitutive equations of the C components for charged species, and the membranepotential (Gawthrop et al., 2017).We believe that thermodynamically consistent models are more likely to remain robustwhen incorporated into whole-cell models and multi-scale models. It is often uncertainwhether models remain predictive under conditions beyond the data they were parameterisedon (Beattie et al., 2018; Fink et al., 2011). As discussed in Smith and Crampin (2004) andTran et al. (2009), many existing models of transporters do not incorporate dependence for allmetabolites and are therefore unable to accurately predict activity under varying metaboliteconcentrations. Because bond graphs enforce physical constraints on models, they provideuseful constraints that dictate the behaviour of the model outside of the range of experimentaldata (Soh and Hatzimanikatis, 2010), being particularly useful if the experimental data werecaptured away from equilibrium. In order to satisfy thermodynamic consistency, bond graphsforce the modeller to incorporate dependence for all metabolites in models of membranetransporters. Adherence to conservation laws may also prove important when transportermodels are incorporated into multi-scale models. In models of cardiac electrophysiology, massand charge conservation impact on long-term behaviour, and violating these conservationlaws may cause the model to drift (Hund et al., 2001; Livshitz and Rudy, 2009; Pan et al.,2018a). Thus in multi-scale models where many individual submodules are simulated overtime-scales far greater than they were originally developed for, accounting for conservationlaws such as conservation of energy may be an important factor in maintaining the stabilityof these models. In this paper, we described the representation of existing models into bond graphs, usinga model of SERCA (Tran et al., 2009) as an example. These methods can be applied toa wide range existing thermodynamically consistent transporter models. A bond graphrepresentation reveals the underlying energetics of these models (Gawthrop and Crampin,2017). In addition to membrane transporters, models of ion channels (Gawthrop et al., 2017),redox reactions (Gawthrop, 2017), metabolic systems (Gawthrop et al., 2015) and signallingpathways (Gawthrop and Crampin, 2016) can also be represented as bond graphs. As thebond graph methodology continues to develop, we anticipate that more biological processessuch as cross-bridge cycling, diffusion and gene regulation will be able to be modelled usingbond graphs. Using the modularity inherent in bond graphs (Gawthrop, 2017), these modelscan be coupled together to assemble whole-cell models that can assess how cells deal withfundamental trade-offs between energy efficiency, speed and robustness.While the bond graph approach can detect thermodynamic inconsistency, it is unable todetect if its thermodynamic parameters are correct. Therefore, in cases such as the model ofNa + /K + ATPase by Terkildsen et al. (2007), parameterising the model with incorrect speciesconstants or equilibrium constants can lead to incorrect equilibrium points even if the modelitself is internally consistent. However, we note that the bond graph framework will flagthese issues to the modeller when two models with conflicting thermodynamic parametersare coupled together. 24 .3. Future work
To date, a large proportion of bond graph models of biochemical systems have beenderived from existing kinetic models (Gawthrop et al., 2015; Gawthrop, 2017; Gawthropand Crampin, 2017). However, as outlined in this paper, developing new models usingthe bond graph framework is a more powerful approach as it guarantees thermodynamicconsistency. A challenge in this process is estimating bond graph parameters from data,as kinetic information alone is insufficient to uniquely determine bond graph parameters(Gawthrop et al., 2015). Therefore future work will focus on using existing information andoptimisation procedures to derive appropriate estimates of bond graph parameters from datawhen thermodynamic data is unavailable, and designing efficient experimental protocols foridentifying bond graph parameters when thermodynamic data is available. In this process,there is a key role for quantifying the uncertainty of biological parameters (Babtie and Stumpf,2017; Beattie et al., 2018). Given that bond graph parameters are physical quantities, itwould be interesting to investigate whether these uncertainties propagate as individual modelsare coupled together.We modelled transporters using multiple states with mass action equations, giving riseto steady-state velocities. However, many models of transporters are simplifications ofthis approach, as they are described by a single equation for the transport rate ratherthan multiple equations for transitions between enzyme states (Keener and Sneyd, 2009).Smith and Crampin (2004) describe a method for reducing multi-state models of membranetransporters into a single equation for transport rate by using rapid equilibrium and steady-state simplifications, while using detailed balance to ensure thermodynamic consistency. Aswe move towards increasingly complex multi-scale models, it is important to develop methodsof simplifying models to reduce computational cost (Smith et al., 2007). An advantage ofusing bond graphs in the simplification process is that they impose constraints to ensurethat simplified models maintain thermodynamic consistency (Gawthrop and Crampin, 2014).
5. Conclusion
We have shown that bond graphs can be used to model membrane transporters whilecapturing fundamental thermodynamic and physical constraints. We apply this frameworkto SERCA and Na + /K + ATPase to develop models that are thermodynamically consistentwhile revealing the underlying energetics of these transporters. Combined with their inherentmodularity, we believe that bond graphs are a powerful tool for incorporating models ofmembrane transporters into other models while maintaining thermodynamic consistency.
Data availability
The code associated with this paper is available from GitHub ( https://github.com/uomsystemsbiology/transporter_thermodynamics ), and is archived on Zenodo ( https://doi.org/10.5281/zenodo.1287353 ) (Pan et al., 2018b). The repositories contain MATLABcode that generates the figures and CellML code containing equations for the models.25 cknowledgements
This research was in part conducted and funded by the Australian Research CouncilCentre of Excellence in Convergent Bio-Nano Science and Technology (project numberCE140100036), and the Australian Research Council’s Discovery Projects funding scheme(project DP170101358). M.P. would like to acknowledge financial support provided by anAustralian Government Research Training Program Scholarship. P.J.G. would like to thankthe Melbourne School of Engineering for its support via a Professorial Fellowship. K.T. issupported by the Heart Foundation of New Zealand (Research Fellowship 1692) and theMarsden Fund Council from Government funding, managed by Royal Society Te Ap¯arangi(Marsden Fast-Start UOA1703).
References
Anton, H., Rorres, C., 2014. Elementary linear algebra : applications version. Hoboken, NJ : John Wiley &Sons Inc.Atkins, P. W., De Paula, J., 2006. Physical chemistry for the life sciences. Oxford University Press ; W.H.Freeman, Oxford, UK : New York.Babtie, A. C., Stumpf, M. P. H., Aug. 2017. How to deal with parameters for whole-cell modelling. Journalof The Royal Society Interface 14 (133), 20170237. doi:10.1098/rsif.2017.0237.Beard, D. A., Qian, H., 2008. Chemical biophysics quantitative analysis of cellular systems. CambridgeUniversity Press, Cambridge; New York.Beattie, K. A., Hill, A. P., Bardenet, R., Cui, Y., Vandenberg, J. I., Gavaghan, D. J., Boer, T. P., Mirams,G. R., Apr. 2018. Sinusoidal voltage protocols for rapid characterisation of ion channel kinetics. TheJournal of Physiology 0 (0). doi:10.1113/JP275733.Blaustein, M. P., Lederer, W. J., Jul. 1999. Sodium/Calcium Exchange: Its Physiological Implications.Physiological Reviews 79 (3), 763–854.Borutzky, W., 2010. Bond Graph Methodology. Springer.de Meis, L., Jul. 2002. Ca -ATPases (SERCA): Energy Transduction and Heat Production in TransportATPases. The Journal of Membrane Biology 188 (1), 1–9. doi:10.1007/s00232-001-0171-5.Fink, M., Niederer, S. A., Cherry, E. M., Fenton, F. H., Koivumäki, J. T., Seemann, G., Thul, R., Zhang, H.,Sachse, F. B., Beard, D., Crampin, E. J., Smith, N. P., Jan. 2011. Cardiac cell modelling: Observationsfrom the heart of the cardiac physiome project. Progress in Biophysics and Molecular Biology 104 (1–3),2–21. doi:10.1016/j.pbiomolbio.2010.03.002.Friedrich, T., Bamberg, E., Nagel, G., Nov. 1996. Na + ,K + -ATPase pump currents in giant excised patchesactivated by an ATP concentration jump. Biophysical Journal 71 (5), 2486–2500. doi:10.1016/S0006-3495(96)79442-0.Gawthrop, P., Bevan, G., Apr. 2007. Bond-graph modeling. IEEE Control Systems 27 (2), 24–45.doi:10.1109/MCS.2007.338279.Gawthrop, P., Smith, L., 1996. Metamodelling: for bond graphs and dynamic systems. Prentice Hallinternational series in systems and control engineering. Prentice Hall, London, New York.Gawthrop, P. J., Apr. 2017. Bond Graph Modeling of Chemiosmotic Biomolecular Energy Transduction.IEEE Transactions on NanoBioscience 16 (3), 177–188. doi:10.1109/TNB.2017.2674683.Gawthrop, P. J., Crampin, E. J., Nov. 2014. Energy-based analysis of biochemical cycles using bond graphs.Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 470 (2171),20140459. doi:10.1098/rspa.2014.0459.Gawthrop, P. J., Crampin, E. J., Mar. 2016. Modular bond-graph modelling and analysis of biomolecularsystems. IET Systems Biologydoi:10.1049/iet-syb.2015.0083.Gawthrop, P. J., Crampin, E. J., Jun. 2017. Energy-based analysis of biomolecular pathways. Proc. R. Soc.A 473 (2202), 20160825. doi:10.1098/rspa.2016.0825. awthrop, P. J., Crampin, E. J., 2018. Biomolecular system energetics. In: Proceedings of the 13thInternational Conference on Bond Graph Modeling (ICBGM’18). Society for Computer Simulation,Bordeaux, available at arXiv:1803.09231.Gawthrop, P. J., Cursons, J., Crampin, E. J., Dec. 2015. Hierarchical bond graph modelling of biochemicalnetworks. Proc. R. Soc. A 471 (2184), 20150642. doi:10.1098/rspa.2015.0642.Gawthrop, P. J., Siekmann, I., Kameneva, T., Saha, S., Ibbotson, M. R., Crampin, E. J., 2017. Bond graphmodelling of chemoelectrical energy transduction. IET Systems Biology 11 (5), 127–138. doi:10.1049/iet-syb.2017.0006.Glitsch, H. G., Jan. 2001. Electrophysiology of the Sodium-Potassium-ATPase in Cardiac Cells. PhysiologicalReviews 81 (4), 1791–1826.Glynn, I. M., 2002. A Hundred Years of Sodium Pumping. Annual Review of Physiology 64 (1), 1–18.doi:10.1146/annurev.physiol.64.081501.130716.Guynn, R. W., Veech, R. L., Oct. 1973. The Equilibrium Constants of the Adenosine Triphosphate Hydrolysisand the Adenosine Triphosphate-Citrate Lyase Reactions. Journal of Biological Chemistry 248 (20),6966–6972.Hansen, P. S., Buhagiar, K. A., Kong, B. Y., Clarke, R. J., Gray, D. F., Rasmussen, H. H., Nov.2002. Dependence of Na + -K + pump current-voltage relationship on intracellular Na + , K + , and Cs + in rabbit cardiac myocytes. American Journal of Physiology - Cell Physiology 283 (5), C1511–C1521.doi:10.1152/ajpcell.01343.2000.Hund, T. J., Kucera, J. P., Otani, N. F., Rudy, Y., Dec. 2001. Ionic Charge Conservation and Long-Term Steady State in the Luo–Rudy Dynamic Cell Model. Biophysical Journal 81 (6), 3324–3331.doi:10.1016/S0006-3495(01)75965-6.Hung, C.-h., Markham, T. L., 1975. The Moore-Penrose inverse of a partitioned matrix. Linear Algebra andits Applications 11 (1), 73–86. doi:10.1016/0024-3795(75)90118-4.Hwang, S.-T., Apr. 2004. Nonequilibrium thermodynamics of membrane transport. AIChE Journal 50 (4),862–870. doi:10.1002/aic.10082.Kawase, Y., Hajjar, R. J., 2008. The cardiac sarcoplasmic/endoplasmic reticulum calcium ATPase: apotent target for cardiovascular diseases. Nature Clinical Practice Cardiovascular Medicine 5 (9), 554–565.doi:10.1038/ncpcardio1301.Keener, J., Sneyd, J., 2009. Mathematical Physiology. Vol. 8/1 of Interdisciplinary Applied Mathematics.Springer New York, New York, NY.Kennedy, J., Eberhart, R., Nov. 1995. Particle swarm optimization. In: IEEE International Conference onNeural Networks, 1995. Proceedings. Vol. 4. pp. 1942–1948 vol.4.Liebermeister, W., Uhlendorf, J., Klipp, E., Jun. 2010. Modular rate laws for enzymatic re-actions: thermodynamics, elasticities and implementation. Bioinformatics 26 (12), 1528–1534.doi:10.1093/bioinformatics/btq141.Livshitz, L., Rudy, Y., Sep. 2009. Uniqueness and Stability of Action Potential Models during Rest,Pacing, and Conduction Using Problem-Solving Environment. Biophysical Journal 97 (5), 1265–1276.doi:10.1016/j.bpj.2009.05.062.Luo, C. H., Rudy, Y., Jun. 1994. A dynamic model of the cardiac ventricular action potential. I.Simulations of ionic currents and concentration changes. Circulation Research 74 (6), 1071–1096.doi:10.1161/01.RES.74.6.1071.Makinose, M., Hasselbach, W., 1971. ATP synthesis by the reverse of the sarcoplasmic calcium pump. FEBSLetters 12 (5), 271–272. doi:10.1016/0014-5793(71)80196-5.Mueckler, M., Thorens, B., Apr. 2013. The SLC2 (GLUT) family of membrane transporters. MolecularAspects of Medicine 34 (2), 121–138. doi:10.1016/j.mam.2012.07.001.Nakao, M., Gadsby, D. C., Sep. 1989. [Na] and [K] dependence of the Na/K pump current-voltage re-lationship in guinea pig ventricular myocytes. The Journal of General Physiology 94 (3), 539–565.doi:10.1085/jgp.94.3.539.Oster, G. F., Perelson, A. S., Katchalsky, A., Feb. 1973. Network thermodynamics: dynamic modelling ofbiophysical systems. Quarterly Reviews of Biophysics 6 (01), 1–134. doi:10.1017/S0033583500000081. an, M., Gawthrop, P. J., Tran, K., Cursons, J., Crampin, E. J., 2018a. Bond graph modelling of thecardiac action potential: Implications for drift and non-unique steady states. Proc. R. Soc. A(In press).doi:10.1098/rspa.2018.0106.Pan, M., Gawthrop, P. J., Tran, K., Cursons, J., Crampin, E. J., Jun. 2018b. Supporting code for “Athermodynamic framework for modelling membrane transporters”.URL https://zenodo.org/record/1287353 Paynter, H. M., 1961. Analysis and design of engineering systems. MIT press.Periasamy, M., Kalyanasundaram, A., Apr. 2007. SERCA pump isoforms: Their role in calcium transportand disease. Muscle & Nerve 35 (4), 430–442. doi:10.1002/mus.20745.Pinz, I., Tian, R., Belke, D., Swanson, E., Dillmann, W., Ingwall, J. S., Mar. 2011. Compromised MyocardialEnergetics in Hypertrophied Mouse Hearts Diminish the Beneficial Effect of Overexpressing SERCA2a.Journal of Biological Chemistry 286 (12), 10163–10168. doi:10.1074/jbc.M110.210757.Polettini, M., Esposito, M., Jul. 2014. Irreversible thermodynamics of open chemical networks. I. Emergent cy-cles and broken conservation laws. The Journal of Chemical Physics 141 (2), 024117. doi:10.1063/1.4886396.Schramm, M., Klieber, H. G., Daut, J., Dec. 1994. The energy expenditure of actomyosin-ATPase, Ca -ATPase and Na + ,K + -ATPase in guinea-pig cardiac ventricular muscle. The Journal of Physiology 481 (Pt3), 647–662.Smith, N. P., Crampin, E. J., Jun. 2004. Development of models of active ion transport for whole-cellmodelling: cardiac sodium–potassium pump as a case study. Progress in Biophysics and Molecular Biology85 (2–3), 387–405. doi:10.1016/j.pbiomolbio.2004.01.010.Smith, N. P., Crampin, E. J., Niederer, S. A., Bassingthwaighte, J. B., Beard, D. A., May 2007. Computationalbiology of cardiac myocytes: proposed standards for the physiome. Journal of Experimental Biology210 (9), 1576–1583. doi:10.1242/jeb.000133.Soh, K. C., Hatzimanikatis, V., Jun. 2010. Network thermodynamics in the post-genomic era. CurrentOpinion in Microbiology 13 (3), 350–357. doi:10.1016/j.mib.2010.03.001.Terkildsen, J., 2006. Modelling Extracellular Potassium Accumulation in Cardiac Ischaemia. Master’s thesis,The University of Auckland.Terkildsen, J. R., Crampin, E. J., Smith, N. P., Nov. 2007. The balance between inactivation and activa-tion of the Na + -K + pump underlies the triphasic accumulation of extracellular K + during myocardialischemia. American Journal of Physiology - Heart and Circulatory Physiology 293 (5), H3036–H3045.doi:10.1152/ajpheart.00771.2007.Thoma, J., Bouamama, B. O., 2000. Modelling and Simulation in Thermal and Chemical Engineering.Springer Berlin Heidelberg, Berlin, Heidelberg.Tran, K., Loiselle, D. S., Crampin, E. J., Jul. 2015. Regulation of cardiac cellular bioenergetics: mechanismsand consequences. Physiological Reports 3 (7), e12464. doi:10.14814/phy2.12464.Tran, K., Smith, N. P., Loiselle, D. S., Crampin, E. J., Mar. 2009. A Thermodynamic Model of theCardiac Sarcoplasmic/Endoplasmic Ca (SERCA) Pump. Biophysical Journal 96 (5), 2029–2042.doi:10.1016/j.bpj.2008.11.045.Walker, M., Williams, G. B., Kohl, T., Lehnart, S., Jafri, M. S., Greenstein, J., Lederer, W. J., Winslow,R., Dec. 2014. Superresolution Modeling of Calcium Release in the Heart. Biophysical Journal 107 (12),3018–3029. doi:10.1016/j.bpj.2014.11.003.Williams, G. B., Chikando, A., Tuan, H.-T. M., Sobie, E., Lederer, W. J., Jafri, M. S., Sep. 2011. Dy-namics of Calcium Sparks and Calcium Leak in the Heart. Biophysical Journal 101 (6), 1287–1296.doi:10.1016/j.bpj.2011.07.021. ppendix A. TF components in biochemical networks For some reactions such as dimerisation, two or more of the same species may be involvedin one side of the reaction. To represent this, the bond connected to that species is fedthrough a TF component. The TF component represents a transformer, which can transmitand convert energy while maintaining energy conservation. The equation for the transformeris e = ne (A.1) f = nf (A.2)Where n represents the stoichiometry of the species. Since e f = e f , the transformer isenergy-conserving as the power in is equal to the power out of the transformer.When applied the biochemical reaction 2 A (cid:10) B , with the bond graph representation inFigure A.1, the constitutive equations of the Re and C components are µ A = RT ln( K A x A ) (A.3) µ B = RT ln( K B x B ) (A.4) v = κ ( e A f /RT − e A r /RT ) (A.5)The equation corresponding to the transformer (with n = 2) is A f = 2 µ A (A.6)˙ x A = − v (A.7)thus the transformer states that the rate of consumption of A is twice the rate of reaction.Because two molecules of A are consumed by the reaction, the forward affinity is twice thechemical potential of A. By combining these equations the law of mass action can be derived: v = κ ( e µ A /RT − e µ B /RT ) = κK A x A − κK B x B (A.8)˙ x A = − v = 2 κK B x B − κ K A x A (A.9)˙ x B = v = κK A x A − κK B x B (A.10) C:A C:BRe:r1TF:n
Figure A.1: Bond graph of the reaction A (cid:10) B . ppendix B. Converting existing kinetic models to bond graphs Because there are a number of existing models of membrane transporters that arethermodynamically consistent, the ability to convert kinetic models of transporters into bondgraphs would save a great amount of effort from having to build new bond graph modelsto explain existing data. In § 2.1, we outlined how a chemical reaction network can berepresented as a bond graph. While the structure of the system can be naturally translateda bond graph, the parameters of existing models do not directly translate into bond graphparameters K and κ . While bond graph parameters define a unique set of kinetic parameters,identifying bond graph parameters from a set of kinetic parameters is more involved. In thissection, we outline how to infer bond graph parameters from kinetic parameters, and discussthe issues that arise during this process.We first use the enzyme cycle example in § 2.1 to illustrate the essential concepts ofinferring bond graph parameters. The relationship between the kinetic and bond graphparameters is given by Eqs. 17–20. Additional constraints are enforced by the detailedbalance constraint in Eq. 33. Therefore, given a set of kinetic constants k +1 , k +2 , k − , k − , anequivalent set of bond graph parameters can be derived by solving the Eqs. 17–20, 33. If wedefine Ln as the element-wise logarithm operator, these constraints can be packaged intothe matrix equation Ln k +1 k +2 k − k − = − Ln κ κ K Si K Se K E K E (B.1)which can be solved through linear methods (Anton and Rorres, 2014). However, it should benoted that the rows of the right hand side are linearly dependent, as R R − R − R − R h − − − i Ln h k +1 k +2 k − k − i T = 0 (B.2)which is easily seen to be the detailed balance constraint in Eq. 34. As a result, a setof bond graph parameters exist for this kinetic system only if the kinetic parameters arethermodynamically consistent, reflecting the fact that bond graphs can only represent thesubset of kinetic systems that are thermodynamically consistent.The equations shown above can be generalised to all chemical reaction networks describedby mass action. The relationships between kinetic and bond graph parameters are capturedin the general equation (Gawthrop et al., 2015) Ln ( k ) = MLn ( W λ ) (B.3)30ith partitions defined as k = k + k − K c , M = I n r × n r N f T I n r × n r N rT N cT , λ = " κK (B.4)Here, k + and k − are vectors of the forward and reverse kinetic rate constants respectively, I n × n is an identity matrix of length n , n r is the number of reactions, N f and N r are theforward and reverse stoichiometric matrices respectively, κ is a vector of the reaction rateconstants, K is a vector of the species thermodynamic constants, K c is a vector of knownequilibrium constants between the species defined in the matrix N c , and W is a diagonalmatrix that accounts for differences in volume between compartments. Therefore the firsttwo rows represent the relationships between kinetic and bond graph constants, and the finalrow represents additional information that is known about the thermodynamic constants.For the enzyme cycle, Eq. B.3 gives rise to Eq. B.1 through the substitutions k + = " k +1 k +2 , k − = " k − k − , κ = " κ κ , K = K Si K Se K E K E , K c = 1 (B.5) N f = , N r = , N c = − , W = I × (B.6)As seen in the enzyme cycle example, Eq. B.3 is not always solvable because there mayexist sets of kinetic rate constants that do not satisfy detailed balance constraints, and aretherefore thermodynamically inconsistent. Because the bond graph approach forces themodeller to use extra discipline in accounting for energy transfer, the exercise of convertingkinetic parameters into bond graphs can be used to verify the thermodynamic consistency ofa kinetic model. We show here that bond graph parameters can be derived from a set ofkinetic parameters if they satisfy a thermodynamic constraint, known in the literature asWegscheider conditions (Liebermeister et al., 2010). By dividing the first n r rows of Eq. B.3by the next n r rows, we find that Ln " K eq K c = " − N T N cT Ln ( K ) (B.7)where N = N r − N f is the stoichiometric matrix. If we define / as the element wise quotientoperator for vectors, K eq = k + /k − is a vector containing the equilibrium constants of thereactions. Thus if Z is a right nullspace matrix (Anton and Rorres, 2014) of [ − N N c ],31ultiplying both sides of Eq. B.7 by Z T on the left gives Z T Ln " K eq K c = 0 (B.8)Thus for any biochemical system with mass action equations, Eq. B.8 can be used to checkthat the kinetic parameters are thermodynamically consistent. As discussed in Gawthropet al. (2015), Eq. B.8 is a form of the Wegscheider condition, which has been associated withthermodynamic consistency (Polettini and Esposito, 2014), therefore only thermodynamicallyconsistent kinetic parameters can be converted into bond graph parameters.Assuming that a solution exists, one solution to Eq. B.3 can be found using the equation λ = W − Exp ( M † Ln ( k )) (B.9)where M † is the Moore-Penrose pseudoinverse of M (Hung and Markham, 1975). If a solutionexists, the set of bond graph parameters is generally not unique. However, as demonstratedin Gawthrop et al. (2015), the affinity and velocity each reaction remains the same regardlessof the parameters chosen.To account for differences in volume between compartments, we include a matrix W inEq. B.3. If a species exists within a compartment, the diagonal entry corresponding to thatspecies is equal to the volume of that compartment. For reaction rate constants and speciesthat do not exist within a volume, the value of the corresponding diagonal element is 1.In Section 3, we use the constraints W i K Nai W e K Nae = 1 (B.10) W i K Ki W e K Ke = 1 (B.11) W i K Cai W sr K Casr = 1 (B.12) W i K MgATP W i K MgADP W i K Pi W isr K H = exp − ∆ G RT ! M = 9881 mM (B.13)where ∆ G = 11 . / mol is the standard free energy of MgATP hydrolysis at 310 K(Guynn and Veech, 1973; Tran et al., 2009). We used the volumes W i = 38 . W e =5 .
182 pL, and W sr = 2 .
28 pL (Luo and Rudy, 1994). For the ATP hydrolysis reaction, weused the intracellular volume for MgATP, MgADP and Pi. In the Terkildsen et al. (2007)model, a volume of W isr = W i was used for H + . The Tran et al. (2009) model describes thecountertransport of Ca and H + , however the steady-state equations were derived assumingthat the cytsol and SR shared the same H + concentration. Therefore we chose to use thevolume W isr = W i + W sr .In examples where charged species are involved (such as in § 2.3), the electrical dependencecan be ignored for the purpose of identifying bond graph parameters, since the components32orresponding to the electrical dependence (the linear C component, and TF components)have a natural bond graph representation.The Tran et al. (2009) and Terkildsen et al. (2007) models assumed that some reactionswere at rapid equilibrium, therefore we approximated them in the bond graph model byreplacing the equilibrium constants with two sufficiently fast kinetic rate constants with thesame equilibrium constant. Appendix C. Modifications to the Terkildsen et al. model of the Na + /K + AT-Pase
Appendix C.1. Updates to equations
The Terkildsen et al. model is a 15-state model of the Na + /K + ATPase (Figure 6).The model was reduced to a 4-state model using rapid equilibrium (Smith and Crampin,2004), and then further simplified using steady state approximations. We found kinetic andthermodynamic issues in the implementation of this model, and resolved these issues asfollows:1. The equilibrium constants of the chemical reaction network were inconsistent with thenumber of binding sites. If we assume that the binding/unbinding events are identical,the kinetic rate constants are proportional to the number of binding sites (Keener andSneyd, 2009). We have corrected the equilibrium constants, shown in red in Figure 6A.2. The original model by Terkildsen et al. (2007) used an incorrect standard free energyof ∆ G = − . / mol for the hydrolysis of ATP, which resulted in an incorrectequilibrium constant. The authors appeared to adjust ∆ G to a physiological pHrather than a pH of 0. Because ∆ G was used in the detailed balance constraint k +1 k +2 k +3 k +4 K d, Na e ( K d, Na e ) ( K d, K i ) k − k − k − k − K d, Na i ( K d, Na i ) ( K d, K e ) K d, MgATP = exp − ∆ G RT ! , (C.1)the model was thermodynamically inconsistent. This error would cause over a 10 –foldincrease in the equilibrium constant at 310 K, therefore we corrected the value to∆ G = 11 . / mol (Guynn and Veech, 1973; Tran et al., 2009).3. Terkildsen et al. (2007) applied a rapid equilibrium approximation to reduce the 15-statemodel into a 4-state model with modified kinetic constants that were functions ofmetabolite concentrations. However, due to algebraic errors, the equations for someof these modified kinetic rate constants ( α +1 , α +3 , α − and α − ) were incorrect. In ourupdated model, we have corrected the following modified rate constants to the equations33elow: α +1 = k +1 ˜Na i, ˜Na i, ˜Na i, ˜Na i, + (1 + ˜Na i, ) + (1 + ˜K i ) − α +3 = k +3 ˜K e ˜Na e, ˜Na e, + (1 + ˜Na e, ) + (1 + ˜K e ) − α − = k − ˜Na e, ˜Na e, ˜Na e, ˜Na e, + (1 + ˜Na e, ) + (1 + ˜K e ) − α − = k − ˜K i ˜Na i, ˜Na i, + (1 + ˜Na i, ) + (1 + ˜K i ) − i, = [Na + ] i K d, Na i e ∆ F V m /RT ˜Na i, = [Na + ] i K d, Na i (C.6)˜Na e, = [Na + ] e K d, Na e e (1+∆) zF V m /RT ˜Na e, = [Na + ] e K d, Na e (C.7)˜K i = [K + ] i K d, K i ˜K e = [K + ] e K d, K e (C.8)and ∆ is the charge translocated by reaction R5. α +1 was one of the modified rate constants that was updated, described in Terkildsenet al. (2007) by the incorrect equation α +1 = k +1 ˜Na i, ˜Na i, (1 + ˜Na i, )(1 + ˜Na i, ) + (1 + ˜K i ) − k +1 scaled by the proportion of state 6 with respectto the total amount of states 1 to 6, a correct expression for the modified rate constantcan be derived as follows: α +1 = k +1 x x + x + x + x + x + x = k +1
11 + x /x + x /x + x /x + x /x + x /x = k +1 − i, + 2 ˜Na − i, ˜Na − i, + ˜Na − i, ˜Na − i, + 2 ˜Na − i, ˜Na − i, ˜K i + ˜Na − i, ˜Na − i, ˜K i = k +1 ˜Na i, ˜Na i, ˜Na i, ˜Na i, + (1 + ˜Na i, ) + (1 + ˜K i ) − Appendix C.2. Reparameterising the Terkildsen et al. model
Using the updated equations of the Terkildsen et al. model, and setting ∆ G =11 . / mol for the detailed balance constraint, we reparameterised the model to the samedata that the Terkildsen et al. model was trained on. We minimised an objective functionthat summarised differences between model predictions and data, using similar methods tothose in Terkildsen (2006), with some minor changes:1. For the extracellular K + data of Nakao and Gadsby (1989) (Figure 8B), the weightingfor extracellular potassium above 5.4 mM was increased from 6 × to 15 × to obtain areasonable fit at physiological concentrations.2. To ensure that the magnitude of cycling velocity matched that of Nakao and Gadsby(1989), we chose to fit the curve for [Na + ] e = 150 mM (Figure 7A) without normalisation.3. Instead of using a local optimiser to minimise the objective function, we used particleswarm optimisation (Kennedy and Eberhart, 1995) followed by a local optimiser tofind a global minimum. Appendix D. Stoichiometry and parameters for SERCA
Table D.1: Forward stoichiometric matrix for the Tran et al. (2009) SERCA model. R1 → → →
2a R4 → → → → →
10 R10 → + able D.2: Reverse stoichiometric matrix for the Tran et al. (2009) SERCA model. R1 → → →
2a R4 → → → → →
10 R10 → + Table D.3: Parameters for the Tran SERCA model.
Adapted from Tran et al. (2009).
Parameter Value Modified kinetic constants k +1 − s − - k +2 − - k +3 . − - k − − - k − − s − - k +3 − s − - K eqd,Cai k + = 3 . × mM − s − k − = 2 . × s − K eqd,Casr k + = 1 . × s − k − = 2 . × mM − s − K eqd,H1 . × − mM k + = 2 . × mM − s − k − = 2 . × s − K eqd,Hi . × − mM k + = 2 . × s − k − = 7 . × mM − s − K eqd,Hsr . × − mM k + = 2 . × mM − s − k − = 2 . × s − K eqd,H . × − mM k + = 2 . × s − k − = 3 . × mM − s − n able D.4: Parameters for the bond graph model of the SERCA pump. The reaction parametershave unit fmol / s, and the species parameters have unit fmol − . Species/Reaction Parameter Value R1 → κ . → κ . → a κ . → κ . → κ . → κ . → κ . → κ . → κ . K . K . K . K . K . K . K . K . K . + K H . i K Cai . sr K Casr . K MgATP . K MgADP . × − P i K Pi . ppendix E. Stoichiometry and parameters for Na + /K + ATPase
Table E.1: Forward stoichiometric matrix for the Terkildsen et al. (2007) model of Na + /K + ATPase.
Reactions1 2 3 4 5 6 7 8 9 10 11 12 13 14 15P1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0P2 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0P3 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0P4 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0P5 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0P6 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0P7 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0P8 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0P9 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0P10 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0P11 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0P12 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0P13 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0P14 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0P15 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1K +i +e +i +e + able E.2: Reverse stoichiometric matrix for the Terkildsen et al. (2007) model of Na + /K + ATPase.
Reactions1 2 3 4 5 6 7 8 9 10 11 12 13 14 15P1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1P2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0P3 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0P4 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0P5 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0P6 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0P7 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0P8 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0P9 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0P10 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0P11 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0P12 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0P13 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0P14 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0P15 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0K +i +e +i +e + able E.3: Kinetic parameters for the updated cardiac Na + /K + ATPase model.
Parameter Description Value k +1 Forward rate constant of reaction R6 1423 . − k − Reverse rate constant of reaction R6 225 . − k +2 Forward rate constant of reaction R7 11564 . − k − Reverse rate constant of reaction R7 36355 . − k +3 Forward rate constant of reaction R13 194 . − k − Reverse rate constant of reaction R13 281037 . − s − k +4 Forward rate constant of reaction R15 30629 . − k − Reverse rate constant of reaction R15 1 . × s − K Voltage-dependent dissociation constant ofintracellular Na + . K Voltage-dependent dissociation constant ofextracellular Na + . K d,Nai Voltage-independent dissociation constant ofintracellular Na + . K d,Nae Voltage-independent dissociation constant ofextracellular Na + . K d,Ki Dissociation constant of intracellular K + .
976 mM K d,Ke Dissociation constant of extracellular K + . K d,MgATP Dissociation constant of MgATP 140 . − . µ m . µ m − able E.4: Parameters for the bond graph version of the updated cardiac Na + /K + ATPasemodel.
Component Description Parameter Value
R1 Reaction R1 κ . / sR2 Reaction R2 κ . / sR3 Reaction R3 κ . / sR4 Reaction R4 κ . / sR5 Reaction R5 κ . / sR6 Reaction R6 κ . / sR7 Reaction R7 κ . / sR8 Reaction R8 κ . / sR9 Reaction R9 κ . / sR10 Reaction R10 κ . / sR11 Reaction R11 κ . / sR12 Reaction R12 κ . / sR13 Reaction R13 κ . / sR14 Reaction R14 κ . / sR15 Reaction R15 κ . / sP Pump state ATP–E i K K . − P Pump state ATP–E i K K . − P Pump state ATP–E i K . − P Pump state ATP–E i Na K . − P Pump state ATP–E i Na K . − P Pump state ATP–E i Na K . − P Pump state P–E i (Na ) K . − P Pump state P–E e Na K . − P Pump state P–E e Na K . − P Pump state P–E e Na K . − P Pump state P–E e K . − P Pump state P–E e K K . − P Pump state P–E e K K . − P Pump state E e (K ) K . − P Pump state ATP–Ee (K ) K . − Ki Intracellular K +i K Ki . − Ke Extracellular K +e K Ke . − Nai Intracellular Na +i K Nai . − Nae Extracellular Na +e K Nae . − MgATP Intracellular MgATP K MgATP . − MgADP Intracellular MgADP K MgADP . × − fmol − P i Free inorganic phosphate K Pi . − H Intracellular H + K H . − mem Membrane capacitance C m z − . z − .9450