BondGraphs.jl: composable energy-based modelling in systems biology

Abstract Summary BondGraphs.jl is a Julia implementation of bond graphs. Bond graphs provide a modelling framework that describes energy flow through a physical system and by construction enforce thermodynamic constraints. The framework is widely used in engineering and has recently been shown to be a powerful approach for modelling biology. Models are mutable, hierarchical, multiscale, and multiphysics, and BondGraphs.jl is compatible with the Julia modelling ecosystem. Availability and implementation BondGraphs.jl is freely available under the MIT license. Source code and documentation can be found at https://github.com/jedforrest/BondGraphs.jl.

Bond graphs represent all systems and reactions in terms of energy flow between elements of a system.Physical variables are described in terms of abstracted efforts (forces, voltages, chemical potentials) and flows (velocities, currents, molar flow rate).Efforts are always in terms of Joules per some quantity x, and flows are always x per second, so that J/x • x/s = J/s, or energy change over time.By abstracting these units, any of the above efforts and flows can be combined in a single model that conserves energy across all domains.
Bond graphs are built with collections of components.Components are generalised forms of physical laws and typically represent a physical element of the system, such as a resistor, point mass, or chemical species.Each component stores a constitutive relation that describes the relationship between the energetic effort [Joules/x] and energy flow [x/second] in the component, where x is a measured quantity (displacement, charge, moles).BondGraphs.jl includes a library of common bond graph component types.These include the C type components (capacitors, chemical species), R type components (resistors, reactions), I type components (inductors, point masses), and energy sources (batteries, chemostats).For more examples and details of bond graph components, refer to Gawthrop and Bevan (2007) and Borutzky (2011).
The biochemical reaction component Re differs from the standard R component in that it is a non-linear function of two efforts (chemical potentials).In other words, a distinction is made between the potential of the reactants and the potential of the products.In bond graph terms, the Re component has two ports instead of one.In other aspects the analogy with the resistive R component still holds, as they are all energy 'sinks' through the release of thermal energy.
The standard constitutive equations for a biochemical bond graph produce mass-action rate laws.Other rate laws may be specified.In the BondGraphs.jldocumentation there is an example of implementing a Michaelis-Menten rate law in place of mass-action kinetics (https://jedforrest.github.io/BondGraphs.jl/stable/examples/).Any rate law can be used in principle, so long as it is thermodynamically consistent.
Components are joined with a bond, representing the energy transfer between components.Each bond has an associated effort and flow variable and a direction that determines sign convention.Bonds may also connect to junctions, which describe particular conservation laws.These are either Equal Effort (effort is equal in connected components) or Equal Flow (flow is identical).In bond graph notation, these are equivalent to the 0 and 1 junctions respectively.Intuitively, these can be thought of as junctions that conserve the (typically) 0D effort and 1D flow.

Thermodynamic consistency
Since bonds are two-way energy connectors, the energy of the entire system is always conserved.(Assuming a closed system.In an open system, energy can be dissipated or added, but it can always be tracked.)A similar rule applies to mass conservation or charge conservation.Bond graphs therefore by design enforce physical and thermodynamic constraints.More precisely, bond graphs satisfy energy conservation as all bonds and junctions are power conserving.While entropy is not always explicitly shown, the components can be modified to account for temperature and entropy flow (for example in Borutzky, 2011).However, it is up to the modeller to ensure that the components are 'well formed' and physically plausible.In other words, if an arbitrary component's constitutive law does not conserve energy (for example, through the use of irreversible rate laws), then the system as a whole will not conserve energy.Nonetheless, the bonds, ports, and junctions in the bond graph framework are thermodynamically safe.Thus, model composition across scales and domains will not break physical constraints.This makes it possible to construct realistic large-scale biophysical models.
Therefore, there is some onus on the modeller to ensure that newly designed components are physically plausible.However, existing standard components (such as the Reaction component Re) are physically plausible.Pre-made components in the BondGraphs.jllibrary all satisfy the laws of thermodynamics.
Specific to biochemistry, bond graphs are able to automatically satisfy detailed balance constraints in chemical reaction network models.This is in contrast to typical mass action models with kinetic parameters, where additional detailed balance constraints need to be identified and enforced (Liebermeister and Klipp, 2006).Bond graphs assign a parameter for each species (thermodynamic constant Ki[mol −1 ]) and reaction (reaction rate κj [mol • s −1 ]) in the system.In contrast, kinetic parameters in conventional models are combinations of thermodynamic parameters, and changing the value of one kinetic parameter without accounting for changes in the other kinetic rates may result in a thermodynamically infeasible system.Previous bond graph papers have addressed this issue (Pan et al., 2019;Gawthrop et al., 2021) and Mason and Covert (2019) benchmark the fitting of thermodynamic parameters to data compared to conventional kinetic parameters.
Bond graphs have the additional benefit of reducing the parameter space to thermodynamically feasible regions.In this sense, they have a similar benefit to convenience kinetics and similar frameworks.A discussion about convenience kinetics and bond graphs is given in Pan et al. (2021).
For further reading on how bond graphs satisfy detailed balance, see Pan et al. (2019).For a discussion on dealing with uncertainty in thermodynamic parameters, we refer to Gawthrop et al. (2015).

B. Model descriptions
Model specifications and code for the two models used in Figure 1 in the Section 3 of the main text: the ion pore transport model (Cudmore et al., 2021), and; the Sarco/Endoplasmic Reticulum Ca2 + -ATPase (SERCA) pump (Pan et al., 2019).
We have included an additional example on a large (281 species, 566 reactions) Saccharomyces cerevisiae (yeast) metabolic network (Stanford et al., 2013).This example features automatically generating a bond graph from a published systems biology model stored in BioModels (Malik-Sheriff et al., 2020).
Each section includes the generated ODE system of equations (or a subset of the equations), model parameters and initial conditions (where applicable), and the Julia code used to create the bond graph models and plots in Figure 1.

B.2.3. Parameters and Initial Conditions
Table 3. Parameter values and initial conditions for the states used in the SERCA bond graph model.Taken from Pan et al. (2019).The input function for Casr was chosen to resemble a physiological calcium spike, characterised by an initial spike in the concentration of Ca 2+ followed by an exponential decay.

B.3. Saccharomyces cerevisiae Metabolic Network
This example demonstrates how one could generate a bond graph model of a large biological network stored in a computer-readable format.For this example, we have used the Saccharomyces cerevisiae metabolic network from Stanford et al. (2013).This model contains 281 species and 566 reactions.The SBML model file contains information on relevant chemical species, reaction stoichiometry, default parameter values, and various other metadata.The file is downloadable from https://doi.org/10.1371/journal.pone.0079195.s002.
The bond graph is constructed using the chemical reaction network interface in BondGraphs.jl.The bond graph in total contains 1019 components, 1127 junctions, and 3813 bonds.We assume mass-action rate laws for simplicity and demonstration purposes; in practice one may update the Re reaction components with more sophisticated rate laws.

B.3.1. ODE
The system of ODEs generated by BondGraphs.jl has 293 equations.For the sake of conciseness, only one representative equation is shown here.
Si are the species amounts, Ki are the thermodynamic constants, rj are the reaction rates.Species and reactions are indexed according to Stanford et al. (2013).

Fig. 1 :
Fig. 1: Bond graph model of a yeast metabolic network (taken from Stanford et al., 2013).Left: 2146×2146 adjacency matrix for the yeast network bond graph.Axes correspond to the indices ai,j of the adjacency matrix, where ai,j = 1 (black dot) means node i (component or junction) is connected to node j.Right: Graph plot of the yeast network bond graph.The graph contains 1019 components (blue for species, orange for reactions and transformers), 1127 junctions (white), and 3813 bonds.

Table 1 .
Bond graph variables and analogous components for common physical domains.