Bond Graph Representation of Chemical Reaction Networks
BBond Graph Representation ofChemical Reaction Networks
Peter J. Gawthrop ∗ and Edmund J. Crampin Systems Biology Laboratory, Department of Biomedical Engineering,Melbourne School of Engineering, University of Melbourne, Victoria 3010,Australia. Systems Biology Laboratory, School of Mathematics and Statistics,University of Melbourne University of Melbourne, Victoria 3010October 16, 2018
Abstract
The Bond Graph approach and the Chemical Reaction Network approach tomodelling biomolecular systems developed independently. This paper brings togetherthe two approaches by providing a bond graph interpretation of the chemical reactionnetwork concept of complexes. Both closed and open systems are discussed.The method is illustrated using a simple enzyme-catalysed reaction and a trans-membrane transporter. ∗ Corresponding author. [email protected] a r X i v : . [ q - b i o . M N ] O c t ontents Introduction
The bond graph method for modelling engineering systems [1–6] was shown to provide athermodynamically consistent approach to modelling biomolecular systems by Oster et al.[7, 8] and further developed by Gawthrop and Crampin [9, 10, 11]. In this context, therelationship between biomolecular systems and electrical circuit theory was explored byOster and Perelson [12].In parallel with the seminal work of Oster et al. [7, 8], the mathematical foundationsof chemical reaction networks (CRN) were being laid by Feinberg [13], Horn and Jackson[14] and Feinberg and Horn [15]. This approach to chemical reaction network theory wasfurther developed by Sontag [16], Angeli [17], and van der Schaft et al. [18, 19, 20]. Generalresults on stability of both closed and open systems of chemical reactions have been derivedand applied to reveal dynamic features of complex (bio)chemical networks [21], dissipationin noisy chemical networks Polettini et al. [22], metabolic networks [23] and multistabilityin interferon signalling Otero-Muras et al. [24].As an energy-based method, bond graphs are related to port-Hamiltonians [25–27]. Aport-Hamiltonian interpretation of CRNs has been given by van der Schaft et al. [28] andthis provides another link between CRNs and bond graphs.The formal concept of complexes is essential to chemical reaction network theory. Com-plexes are the combination of chemical species forming the substrate and products of thenetwork reactions. This paper links chemical reaction network theory to the bond graphapproach by incorporating the concept of complexes into bond graph modelling of biomolec-ular systems. § § § § § The notion of complexes was defined by Feinberg and Horn [15]: “By the complexes in amechanism we mean the set of entities appearing before or after arrows in that mechanism.”where “mechanism” is a generalisation of “chemical reaction”. This section introduces somebasic ideas relating to the use of complexes in describing chemical reaction networks bymeans of the simple reaction network exampleA + E r C r B + E (1)This example involves the four species A, B, C and E and the two reactions r and r .It represents the reaction A B catalysed by the enzyme E and with intermediatecomplex C [29]. The substrate of reaction r is A + E and the product is C; the substrateof reaction r is C and the product is B + E. Thus there are three complexes associated3ith this reaction network: A + E, C and B + E; C forms not only the right-hand side ofreaction r but also the left-hand side of reaction r .The standard stoichiometric approach would be to define the species state x and reac-tion flow v as: x = x A x B x C x E v = (cid:18) v v (cid:19) (2)where v and v are the flows associated with reactions r and r respectively. The rate ofchange of species ˙ x is then given in terms of the stoichiometric matrix N and reaction flow v as : ˙ x = N v where N = − − − (3)In contrast, the complex-based approach uses the complex flows v c as an intermediatequantity. Thus define v c = v c v c v c (4)where v c , v c and v c are the flows associated with complexes A + E, C and B + E respec-tively.The rate of change of species ˙ x is given in terms of the matrix Z and the complex flow v c as: ˙ x = Zv c where Z = (5)and the complex flow v c is given in terms of the matrix D and the reaction flow v as: v c = Dv where D = − −
10 1 (6)If follows from Equations (16) and (18) that ˙ x = ZDv and thus it follows from (3) that N = ZD (7) Although equation (3) is linear, as discussed in § v is, in general, a nonlinear function of the species state x . In this particular case the expression for v involves the nonlinear terms x E x A and x E x B . +E Cr1 B+Er2 Figure 1: Digraph corresponding to the D matrix (6) for the system A + E CB + E. The three complexes A + E, C and B + E appear as nodes connected by a digraphwith edges corresponding to the two reactions r and r .The fundamental motivation for the complex-based approach is that graph theory canbe applied to the directed graph formed by taking the complexes to be vertices and thereactions to be edges. In particular, D is the incidence matrix of the graph and has theproperty that each column of D contains exactly one 1 and exactly one −
1; the otherelements being zero. The corresponding digraph (plotted using graphviz ) [30]) appears inFigure 1.Following Gawthrop and Crampin [9], the stoichiometric matrix N can be written as: N = N r − N f (8)where N f and N r connect the forward and reverse sides of the reaction to species . In asimilar fashion, D can be written as: D = D r − D f (9)where D f and D r connect the forward and reverse sides of the reaction to complexes . D f and D r can always be deduced from D as D f and D r correspond to the negative andpositive elements of D respectively.The columns of N f correspond to the substrate complexes and that the columns of N r correspond to the product complexes. It follows that the columns of both of these matricescontain all of the relevant complexes, possibly repeated. Hence Z can be obtained asfollows:1. Create the matrix Z from N f and N r and create the corresponding matrix D Z = (cid:16) N f ... N r (cid:17) (10) D = − I n V × n V . . .I n V × n V (11)It follows from Equation (8) that Z D = N .2. Delete repeated columns of Z to create Z and sum the corresponding rows of D tocreate D . 5ontinuing the example of this section N f = N r = (12)and so Z = D = − −
11 00 1 (13)As columns two and three are identical, column three of Z is deleted to give Z (5), androws two and three of D are merged to give D (6). Figure 2(a) shows the approach used by Gawthrop and Crampin [9, 10, 11] to repre-sent closed systems . However, following the approach of Gawthrop [31], the
Faraday-equivalent potential φ , with units of V, is used in place of chemical potential µ with unitsof J mol − . C represents the n s C components representing the chemical species; φ is thevector of the chemical potentials, and ˙ x the corresponding flow rates. R e represents the n r Re components representing the chemical reactions with forward and reverse potentialΦ f and Φ r and flow rate v . T F : N r and T F : N r represent the bond graph transformers encapsulating the system stoichiometry. A key feature of transformers is that they relateboth the efforts and flows on the corresponding bonds whilst conserving energy [9–11].Thus, with reference to Figure 2(a)˙ x = ˙ x r − ˙ x f = N r v − N f v = N v (14)Φ = Φ f − Φ r = N f T φ − N rT φ = − N T φ (15)In contrast, Figure 2(b) shows the complex-based approach used here. T F : Z representsthe bond graph transformer relating the n s species to the n c complexes which then becomethe reaction forward complex (substrates) via T F : D f and the reaction reverse complex(products) via T F : D r . With reference to Figure 2(b), the transformer equations become:˙ x = Zv c (16) φ c = Z T φ (17) v c = v cr − v cf = D r v − D f v = Dv (18)Φ = Φ f − Φ r = D f T φ c − D rT φ c = − D T φ c (19)6 φ ˙ xφ ˙ x f φ ˙ x r T F : N r T F : N f CR e v Φ f v Φ r (a) Standard approach O φ c v c φ c v cf φ c v cr T F : D r T F : D f R e v Φ f v Φ r T F : Z C φ ˙ x (b) Complex-based approach Figure 2: A Bond Graph Approach to Complexes. (a) The standard approach given byGawthrop and Crampin [9, 10, 11]. The bond symbols (cid:43) correspond to vectors of bonds; C , R e and O correspond to arrays of C , Re and components; the two T F compo-nents represent the intervening junction structure comprising bonds, and junctionsand TF components. N f and N r are the forward and reverse stoichiometric matrices.(b) The complex based approach. T F : Z represents the junction structure connectingcomplexes and species where Z appears in Equation (16) and T F : D f T F : D r representthe junction structure connecting reactions and complexes where D f and D r appear inEquations (18) and (9). 7 .1 Example: A + E C B + E φ A ˙ x A φ A φ B v v v v v v φ E φ C φ E v φ E v φ B ˙ x B ˙ x c ˙ x E C : A C : B C : E C : C R e : r R e : r (a) Standard approach v c φ c φ c v c φ c v c φ E φ B φ A ˙ x A φ C ˙ x C ˙ x B ˙ x E v v C : B C : E C : A C : C : r1
10 Re : r2
01 00 (b) Complex-based approach
Figure 3: Example: A + E C B + E. (a) A bond graph without explicit repre-sentation of complexes. (b) The complex covariables correspond to the three highlightedbonds. The junction structure connecting the three highlighted bonds to the speciescorresponds to
T F : Z of Figure 2(b) and the junction structure connecting the reaction Re components to the three highlighted bonds corresponds to T F : D f and T F : D r of Fig-ure 2(b). Bonds pointing into the Re components correspond to T F : D f , those pointing away from the Re components correspond to T F : D r .Figure 3(a) gives the standard bond graph for the reaction A + E C B + Ecorresponding to the general Figure 2(a) and Figure 3(b) gives the complex-based bondgraph corresponding to the general Figure 2(b). The three bonds corresponding to thethree complex efforts Φ c and flows v c are highlighted. This section derives the main properties of the CRN modelled by a bond graph includingthe complex concept. The notation and concepts of van der Schaft et al. [20, §
2] are usedand reversible reactions are used at the outset. Following the notation of van der Schaftet al. [20], · , xx (cid:11) , Ln and Exp denote elementwise multiplication, division, natural logarithmand exponentiation of column vectors. In particular for two column vectors x and y : x · y = diag( x ) y xy = (diag( y )) − x (20)The basic equation for the potential of species expressed as the Faraday-equivalent8otential [31] is φ = φ (cid:11) + φ N Ln xx (cid:11) (21)where φ N = RTF ≈
26 mV (22)Alternatively, (21) can be rewritten as φ = φ N Ln ( K s · x ) (23)where K s = exp φ (cid:11) φ N x (cid:11) (24) The basic bond graph notion of transformers as expressed in Figure 2(b) means that thepotential of the complexes can be expressed as:Φ c = Z T φ (25)hence Φ c = Φ c (cid:11) + φ N Ln XX (cid:11) (26)where Φ c (cid:11) = Z T φ (cid:11) (27) X = Exp (cid:0) Z T Ln x (cid:1) = n s (cid:89) j =1 x z ji j (28)and X (cid:11) = Exp (cid:0) Z T Ln x (cid:11) (cid:1) = n s (cid:89) j =1 x (cid:11) jz ji (29)Alternatively, using Equation (23)Φ c = φ N Ln ( K c · X ) (30)where K c = Exp (cid:0) Z T Ln K s (cid:1) = n s (cid:89) j =1 k sj z ji (31)Using Equation (20), Equation (30) can also be written asΦ c = φ N Ln (diag K c X ) (32) Mass action kinetics correspond to the Marcelin-de Donder formula [8, 9, 32]: v = κ · (cid:18) Exp Φ f φ N − Exp Φ r φ N (cid:19) (33)where Φ f = N f T φ = D f T Φ c (34)and Φ r = N rT φ = D rT Φ c (35)9sing Equation (30), (33) becomes: v = κ · (cid:18) Exp D f T Φ c φ N − Exp D rT Φ c φ N (cid:19) (36)Because the matrices D f T and D rT are simply selecting the appropriate complexes for eachreaction, each row has exactly one unit element and the rest zero. Hence Equation (36)becomes: v = κ · (cid:18) D f T Exp Φ c φ N − D rT Exp Φ c φ N (cid:19) = − κ · D T Exp Φ c φ N (37)Using Equation (32) Equation (37) becomes: v = K v X = K v Exp (cid:0) Z T Ln x (cid:1) (38)where K v = − κ · (cid:0) D T diag K c (cid:1) (39)Hence the system state equation for mass action kinetics is:˙ x = ZDK v X = N K v X = N K v Exp (cid:0) Z T Ln x (cid:1) (40)This is essentially Equation (4) of van der Schaft et al. [20]. Note that the term Exp (cid:0) Z T Ln x (cid:1) appearing in equations (38) and (40) is, in general, nonlinear. As will be seen in the fol-lowing section, this term leads to products of species states. Substituting the numerical values from the example of § K c = Exp (cid:0) Z T Ln K s (cid:1) = Exp ln K sA ln K sB ln K sC ln K sE = Exp ln K sA + ln K sE ln K sC ln K sB + ln K sE = K sA K sE K sC K sB K sE (41)10imilarly: X = x A x E x C x B x E (42)Substituting the numerical values from the example of § K v = − κ · (cid:0) D T diag K c (cid:1) = − κ · (cid:18) − − (cid:19) diag K c = (cid:18) κ K sA K sE − κ K sC κ K sC − κ K sB K sE (cid:19) (43)Hence, using (38) v = K v X = (cid:18) κ ( K sA K sE x A x E − K sC x c ) κ ( K sC x c − K sB K sE x B x E ) (cid:19) (44) The seminal book “Free energy transduction and biochemical cycle kinetics” of Hill [33]contains an example of a membrane transporter which is discussed in detail by Gawthropand Crampin [11]. The bond graph is given in Figure 4(a) and the bond graph redrawnto expose the complexes is given in Figure 4(b); the ten bonds corresponding to the tencomplex efforts Φ c and flows v c are highlighted.The ten complexes are: E, Mi + E, EM, Li + EM, LEM, LEsM, Lo + EsM, EsM,Mo + Es and Es. They are connected by the seven reactions em, lem, lesm, esm, es, e andslip. The corresponding digraph (plotted using graphviz ) [30]) appears in Figure 6(a). Enzyme-catalysed reactions such as (1), § κ e given in terms of the rate11 e : e R e : l es m C:MoC:LoC:MiC:Li C:EC:LEMC:EM C:LEsMC:EsMC:Es00 00 R e : s li p (a) Standard approach R e : e R e : l es m Re:lem Re:esmRe:em Re:es R e : s li p (b) Complex-based approach Figure 4: Example: Transporter [33]. (a) The bond graph without explicit representationof complexes [11]. (b) The ten complexes correspond to the ten highlighted bonds. Thejunction structure connecting the ten highlighted bonds to the species corresponds to
T F : Z of Figure 2(b) and the junction structure connecting the reaction Re components to theten highlighted bonds corresponds to T F : D f and T F : D r of Figure 2(b).12onstants κ and κ of the reactions r and r as κ e = e ¯ κK c k m + σ v (45)where k m = K c K e (46)¯ κ = κ κ κ + κ (47)and σ v = (cid:40) exp Φ fRT +exp Φ rRT κ = κ exp Φ f RT κ (cid:29) κ (48)where Φ f and Φ r are the overall forward and reverse reaction potentials and e is the totalamount of enzyme both free and bound to C. In particular, in the case of the reactions of(1): σ v = (cid:40) K A x A + K B x B κ = κ K A x A κ (cid:29) κ (49)When dealing with networks of enzyme catalysed reactions such as (1) there are twochoices: either explicitly model the intermediate species C and use two reactions with constant values of κ or use a single reaction approximation without intermediate speciesC and an equivalent rate-constant κ e which is a function of the species states x .However, as discussed by Gunawardena [34], this approximation should be used withcare to avoid violating the fundamental laws of thermodynamics. For example, whenmodelling networks such as the mitogen-activated protein kinase (MAPK) cascade whereenzymes compete and are themselves reaction products, it has been argued [35, § κ is replaced by an expression κ e ( x ) dependent on species states x . There are a number of ways of converting closed systems to open systems whilst retainingthe basic closed system formulation. Horn and Jackson [14] use the concept of a zerocomplex to act as a generalised source and sink of chemical species and this idea is followedup by van der Schaft et al. [20]. Polettini and Esposito [36] use the concept of a chemostat to act as a source and sink of chemical species at fixed concentration and this idea isfollowed up by Gawthrop and Crampin [10]. The chemostat has three interpretations:13. one or more species is fixed to give a constant concentration [37]; this implies thatan appropriate external flow is applied to balance the internal flow of the species.2. an ideal feedback controller is applied to species to be fixed with setpoint as the fixedconcentration and control signal an external flow.3. as a C component with a fixed state.The chemostat approach is used here.As discussed by Gawthrop and Crampin [10], for each species set to be a chemostat,the corresponding row in the stoichiometric matrix N is replaced by a zero vector to formthe chemodynamic stoichiometric matrix N cd . Using the same motivation as that leadingto equation (7), N cd is written as: N cd = Z cd D cd (50)In this case, the closed-system equations (16)– (19) are replaced by˙ x = Z cd v c (51) φ c = Z T φ (52) v c = v cr − v cf = D cd v (53)Φ = Φ f − Φ r = D f T φ c − D rT φ c = − D T φ c (54)Note that it is the flow equations (51) and (53) that are changed; the potential equations(52) and (54) remain the same as those for the closed system (17) and (19). In particular,some complexes associated with Z and D , and thus the potential equations (52) and (54)are not associated with Z cd and D cd , and thus the flow equations (51) and (53). Hencethe digraph associated with D cd does not necessarily contain all of the complex nodesassociated with D . E Cr1r2
Figure 5: Digraph corresponding to the D cd matrix (56) for the system A + E CB + E of § §
3. Compared to Figure 1, setting the species A and B to be chemostatsreduces the number of complexes to two and the digraph is cyclic.14n the case of the system A + E r C r B + E and choosing the two species Aand B to be chemostats, equation (3) is replaced by:˙ x = N cd v where N cd = − − (55)Thus the two chemostats have constant state x A and x B . The decomposition of Equation(50) gives: Z cd = D cd = (cid:18) − − (cid:19) (56)The digraph corresponding to D cd is given in Figure 5; this corresponds to the flow equa-tions (51) and (53). On the other hand, the digraph corresponding to D is given in Figure 1;this corresponds to the potential equations (52) and (54). Thus the cyclic flow associatedwith the digraph of Figure 5 is driven by the potentials associated with the digraph ofFigure 1. The closed system digraph, corresponding to D and the potential equations of the opensystem, is given in Figure 6(a).As discussed by [11], the open system is created by choosing the four species: Li, Lo,Mi and Mo to be chemostats. The flow digraph with incidence matrix D cd of Figure 6(b)has six nodes corresponding to the complexes: E, EM, LEM, LEsM, EsM and Es. Thisdigraph still has the seven connecting reactions listed in § i+EEMem EsMslipLi+EMLEMlem LEsMlesm Lo+EsMesmMo+EsesEsE e (a) Closed EEMemLEMlem EsMslip LEsMlesm esmEsese (b) Open
Figure 6: Digraphs corresponding to the D matrix (18) for the closed and open systemsfor the transporter system. (a) The ten nodes corresponding to the ten complexes areconnected by three disjoint linear graphs. (b) The four chemostats reduce the number ofcomplexes to six and the corresponding six nodes are connected by a cyclic digraph.16 Conclusion
The complex approach to modelling chemical reaction networks as introduced by Feinberg[13], Horn and Jackson [14] and Feinberg and Horn [15] and expanded by van der Schaftet al. [18, 19, 20, 28] has been given a bond graph interpretation thus enabling resultsfrom the complex approach to be applied to the bond graph approach and vice versa . Inparticular, the decomposition of the stoichiometric matrix N into the complex compositionmatrix [20] Z and the complex graph incidence matrix D (where N = ZD ) is given a bondgraph interpretation.The approach is developed for closed systems, but extended to open systems via thepreviously developed notion of chemostats [10, 36]. The corresponding chemodynamicstoichiometric matrix N cd [10] is decomposed into the chemodynamic complex compositionmatrix Z cd and the chemodynamic complex graph incidence matrix D cd (where N cd = Z cd D cd ). The complex graph incidence matrix D determines both the flow and potentialof closed systems, but in open systems the flow is determined by D cd and the potential by D . As, in general D cd (cid:54) = D , the digraph for the flow of open systems is not the same asthe digraph for potentials. In particular, with reference to Figure 6, the flow and potentialdigraphs for open systems may be structurally different.The combination of the explicit energy-compliance feature of the bond graph modellingapproach with the generic results of the graph-theory based chemical reaction networkapproach will, it is hoped, lead to new results and methods for the analysis and synthesisof biomolecular systems. Peter Gawthrop would like to thank the Melbourne School of Engineering for its supportvia a Professorial Fellowship. This research was in part conducted and funded by theAustralian Research Council Centre of Excellence in Convergent Bio-Nano Science andTechnology (project number CE140100036). The authors would like to thank Ivo Siekmannfor alerting them to references [18–20] and the anonymous reviewers for helpful commentson the manuscript.
References [1] H. M. Paynter.
Analysis and design of engineering systems . MIT Press, Cambridge,Mass., 1961. 3[2] F. E. Cellier.
Continuous system modelling . Springer-Verlag, 1991.[3] P. J. Gawthrop and L. P. S. Smith.
Metamodelling: Bond Graphs and DynamicSystems . Prentice Hall, Hemel Hempstead, Herts, England., 1996. ISBN 0-13-489824-9. 174] Peter J Gawthrop and Geraint P Bevan. Bond-graph modeling: A tutorial introduc-tion for control engineers.
IEEE Control Systems Magazine , 27(2):24–45, April 2007.doi:10.1109/MCS.2007.338279.[5] Wolfgang Borutzky.
Bond graph methodology: development and analysis of multidis-ciplinary dynamic system models . Springer, Berlin, 2010. ISBN 978-1-84882-881-0.doi:10.1007/978-1-84882-882-7.[6] Dean C Karnopp, Donald L Margolis, and Ronald C Rosenberg.
System Dynamics:Modeling, Simulation, and Control of Mechatronic Systems . John Wiley & Sons, 5thedition, 2012. ISBN 978-0470889084. 3[7] George Oster, Alan Perelson, and Aharon Katchalsky. Network thermodynamics.
Nature , 234:393–399, December 1971. doi:10.1038/234393a0. 3[8] George F. Oster, Alan S. Perelson, and Aharon Katchalsky. Network thermodynamics:dynamic modelling of biophysical systems.
Quarterly Reviews of Biophysics , 6(01):1–134, 1973. doi:10.1017/S0033583500000081. 3, 9[9] Peter J. Gawthrop and Edmund J. Crampin. Energy-based analysis of biochemicalcycles using bond graphs.
Proceedings of the Royal Society A: Mathematical, Physicaland Engineering Science , 470(2171):1–25, 2014. doi:10.1098/rspa.2014.0459. Availableat arXiv:1406.2447. 3, 5, 6, 7, 9, 11[10] P. J. Gawthrop and E. J. Crampin. Modular bond-graph modelling and analysis ofbiomolecular systems.
IET Systems Biology , 10(5):187–201, October 2016. ISSN 1751-8849. doi:10.1049/iet-syb.2015.0083. Available at arXiv:1511.06482. 3, 6, 7, 13, 14,17[11] Peter J. Gawthrop and Edmund J. Crampin. Energy-based analysis of biomolecularpathways.
Proceedings of the Royal Society of London A: Mathematical, Physical andEngineering Sciences , 473(2202), 2017. ISSN 1364-5021. doi:10.1098/rspa.2016.0825.Available at arXiv:1611.02332. 3, 6, 7, 11, 12, 15[12] G. Oster and A. Perelson. Chemical reaction networks.
Circuits and Sys-tems, IEEE Transactions on , 21(6):709 – 721, November 1974. ISSN 0098-4094.doi:10.1109/TCS.1974.1083946. 3[13] Martin Feinberg. On chemical kinetics of a certain class.
Archive for Rational Me-chanics and Analysis , 46(1):1–41, 1972. ISSN 0003-9527. doi:10.1007/BF00251866. 3,17[14] F. Horn and R. Jackson. General mass action kinetics.
Archive for Rational Mechanicsand Analysis , 47(2):81–116, Jan 1972. ISSN 1432-0673. doi:10.1007/BF00251225. 3,13, 17 1815] Martin Feinberg and Friedrich J.M. Horn. Dynamics of open chemical systems and thealgebraic structure of the underlying reaction network.
Chemical Engineering Science ,29(3):775 – 787, 1974. ISSN 0009-2509. doi:10.1016/0009-2509(74)80195-8. 3, 17[16] E.D. Sontag. Molecular systems biology and control.
European Journal of Control ,11:1–40, 2006. 3[17] David Angeli. A tutorial on chemical reaction network dynamics.
Eu-ropean Journal of Control , 15(34):398 – 406, 2009. ISSN 0947-3580.doi:http://dx.doi.org/10.3166/ejc.15.398-406. 3[18] A. van der Schaft, S. Rao, and B. Jayawardhana. On the mathematical structure ofbalanced chemical reaction networks governed by mass action kinetics.
SIAM Journalon Applied Mathematics , 73(2):953–973, 2013. doi:10.1137/11085431X. 3, 17[19] Arjan van der Schaft, Shodhan Rao, and Bayu Jayawardhana. Complex and detailedbalancing of chemical reaction networks revisited.
Journal of Mathematical Chemistry ,53(6):1445–1458, Jun 2015. ISSN 1572-8897. doi:10.1007/s10910-015-0498-2. 3, 17[20] A. J. van der Schaft, S. Rao, and B. Jayawardhana. A network dynamics approachto chemical reaction networks.
International Journal of Control , 89(4):731–745, 2016.doi:10.1080/00207179.2015.1095353. 3, 8, 10, 13, 17[21] Carsten Conradi, Dietrich Flockerzi, J¨org Raisch, and J¨org Stelling. Sub-network analysis reveals dynamic features of complex (bio) chemical networks.
Proceedings of the National Academy of Sciences , 104(49):19175–19180, 2007.doi:10.1073/pnas.0705731104. 3[22] M. Polettini, A. Wachtel, and M. Esposito. Dissipation in noisy chemical networks:The role of deficiency.
The Journal of Chemical Physics , 143(18):184103, 2015.doi:10.1063/1.4935064. 3[23] Oleksandr Ivanov, Arjan van der Schaft, and Franz J. Weissing. Steady states andstability in metabolic networks without regulation.
Journal of Theoretical Biology ,401:78 – 93, 2016. ISSN 0022-5193. doi:10.1016/j.jtbi.2016.02.031. 3[24] Irene Otero-Muras, Pencho Yordanov, and Joerg Stelling. Chemical reaction networktheory elucidates sources of multistability in interferon signaling.
PLOS Computa-tional Biology , 13(4):1–28, 04 2017. doi:10.1371/journal.pcbi.1005454. 3[25] G. Golo, A.J. van der Schaft, P.C. Breedveld, and B.M. Maschke. Hamiltonian for-mulation of bond graphs. In R. Johansson and A. Rantzer, editors,
Nonlinear andHybrid Systems in Automotive Control , pages 351–372. Springer, London, 2003. 3[26] D. Vink, D. Ballance, and P. Gawthrop. Bond graphs in model matching control.
Mathematical and Computer Modelling of Dynamical Systems , 12(2-3):249 – 261, 2006.doi:10.1080/13873950500068278. 1927] Alejandro Donaire and Sergio Junco. Derivation of input-state-output port-Hamiltonian systems from bond graphs.
Simulation Modelling Practice and Theory , 17(1):137 – 151, 2009. ISSN 1569-190X. doi:10.1016/j.simpat.2008.02.007. Bond GraphModelling. 3[28] A.J. van der Schaft, S. Rao, and B. Jayawardhana. On the network thermodynamicsof mass action chemical reaction networks.
IFAC Proceedings Volumes , 46(14):24 – 29,2013. ISSN 1474-6670. doi:10.3182/20130714-3-FR-4040.00001. 1st IFAC Workshopon Thermodynamic Foundations of Mathematical Systems Theory. 3, 17[29] James P Keener and James Sneyd.
Mathematical Physiology: I: Cellular Physiology ,volume 1. Springer, 2nd edition, 2009. 3[30] Emden R. Gansner and Stephen C. North. An open graph visualization system andits applications to software engineering.
SOFTWARE - PRACTICE AND EXPERI-ENCE , 30(11):1203–1233, 2000. 5, 11[31] P. J. Gawthrop. Bond graph modeling of chemiosmotic biomolecular energy trans-duction.
IEEE Transactions on NanoBioscience , 16(3):177–188, April 2017. ISSN1536-1241. doi:10.1109/TNB.2017.2674683. Available at arXiv:1611.04264. 6, 9[32] Pierre Van Rysselberghe. Reaction rates and affinities.
The Journal of ChemicalPhysics , 29(3):640–642, 1958. doi:10.1063/1.1744552. 9[33] Terrell L Hill.
Free energy transduction and biochemical cycle kinetics . Springer-Verlag,New York, 1989. 11, 12[34] Jeremy Gunawardena. Time-scale separation – Michaelis and Menten’s old idea,still bearing fruit.
FEBS Journal , 281(2):473–488, 2014. ISSN 1742-4658.doi:10.1111/febs.12532. 13[35] Eberhard O. Voit.
A First Course in Systems Biology . Garland Science, New Yorkand London, 2013. 13[36] Matteo Polettini and Massimiliano Esposito. Irreversible thermodynamics of openchemical networks. I. Emergent cycles and broken conservation laws.
The Journal ofChemical Physics , 141(2):024117, 2014. doi:10.1063/1.4886396. 13, 17[37] Peter J. Gawthrop, Joseph Cursons, and Edmund J. Crampin. Hierarchical bondgraph modelling of biochemical networks.