Energy-based advection modelling using bond graphs

Advection, the transport of a substance by the flow of a fluid, is a key process in biological systems. The energy-based bond graph approach to modelling chemical transformation within reaction networks is extended to include transport and thus advection. The approach is illustrated using a simple model of advection via circulating flow and by a simple pharmacokinetic model of anaesthetic gas uptake. This extension provides a physically consistent framework for linking advective flows with the fluxes associated with chemical reactions within the context of physiological systems in general and the human physiome in particular.


Introduction
Advection is the transport of a substance by the flow of a fluid. Modelling the physiome [1][2][3], in particular modelling circulatory systems [4][5][6], requires methods for modelling the molar flows due to advection. Bond graphs have recently been used in physiome modelling due to their ability to model a wide range of biophysical systems [5]. To date, the bond graph approach to modelling biomolecular systems [7][8][9] considers the transformation of molecules via chemical reactions and the corresponding molar flow rates.
However, as pointed out by Cellier [10, p. 397] 'The term "molar flow" has … been used in two quite different contexts. On the one hand, it denotes the physical transport of matter from one point in space and time to another, while on the other hand, it describes the transformation of one chemical species into another during a chemical reaction.' In fact, both types of molar flow are needed in physiological models. For example, the binding and unbinding of oxygen dissolved in blood to haemoglobin can be viewed as a chemical reaction leading to transformation molar flow; however, the advection of haemoglobin though the blood stream transports the haemoglobin from the lungs to other organs and back again leading to advection molar flow. There is therefore a need to bring both transformation and advection under a common modelling framework; this paper shows how the bond graph approach can provide this framework.
Previously, convection bonds [11,12], which carry two effort variables, have been proposed, and advection has been modelled in the context of chemical reactors by Couenne et al. [13]. By contrast, this paper uses the standard bond graph representation with a single effort variable thus allowing standard bond graph notation and software to be used. Although energy-based modelling of distributed (PDE) systems is possible [14], this paper takes the lumped approach which is typically used for modelling the cardiovascular system [15][16][17] and allows the use of well-established fast ODE solvers. This paper shows how advection can be included within a bond graph model of onedimensional fluid flow. The generation of the underlying fluid flow model is not part of this paper; general considerations such as the number of compartments to use to approximate a more detailed model of haemodynamics taking account of flow characteristics such as the Womersley number are available [16,17]. However, the advection approach of this paper can be straightforwardly coupled to pre-existing fluid dynamics models. Section 2 examines hydrochemical transduction-the transduction of energy between fluid flow and advected chemical potential-from a bond graph perspective. Section 3 looks at the advection of chemical species though an orifice and though a pipe and suggests a new bond graph component: RA. Section 4 compares and contrasts advection of chemical species via the RA component and transformation of chemical species via the bond graph Re component. Section 5 shows how circulatory advection (such as the human blood circulation) combines with binding and unbinding of a ligand and enzyme (such as haemoglobin and oxygen). Section 6 gives a numerical example of a pharmacokinetic system, representing the delivery of a gaseous anaesthetic to a human subject, based on the models developed by Mapleson [18][19][20]. Section 7 summarizes the paper and suggests further research.
The bond graph theory used in this paper was developed by Gawthrop and Crampin [9] and introductory material is available [21,22]. The bond graph modelling in this paper is based on the BondGraphTools Python package [23] available at https://pypi.org/project/BondGraphTools/. The code used for the examples in this paper is available at https://github.com/gawthrop/Advection22.

Hydrochemical transduction
Although different physical and chemical domains have different quantities and units, they share the same energy with units of joules (J). This fact is used by the bond graph approach to provide a unified approach to energy transduction between different energy domains using the transformer TF and gyrator GY elements.
This section looks at hydrochemical transduction with unidirectional incompressible hydraulic flow; bidirectional flow is considered in §3. In the incompressible hydraulic domain, the two bond graph covariables are pressure P [Pa] and volumetric flow Q [m 3 s −1 ] [24]. Noting that the unit Pa can be rewritten as J m −3 , the product of the covariables has units J s −1 . In the chemical domain, the two bond graph covariables are chemical potential μ [J mol −1 ] and advective molar flow 1 f [mol s −1 ]. The flow variables are related by where c i [mol m −3 ] is the concentration of the substance in the liquid. Figure 1 shows the bond graph of a two-domain system. The section labelled hydraulic represents the incompressible hydraulic flow alluded to above. The R:r i component represents an orifice though which the fluid flows; the flow Q is typically a nonlinear function of the net pressure across the orifice [5,24]. The section marked chemical represents a single substance being carried by the fluid though the orifice.
The hydraulic and chemical domains are connected by the transformer TF:c i component with modulus c i , which ensures that the chemical flow f and hydraulic flow Q are related by equation (2.1). The bond graph component transformer TF:c transmits energy, but does not dissipate energy. It follows that ð2:2Þ With reference to figure 1, the energy flow, or power P TF , transferred by the transformer TF:c is the product of the flow and effort variables The power dissipated in the R:r i component is Thus the power dissipated in the R:r i component is the sum of the hydraulic dissipation P H and the power P TF transmitted from the chemical domain.

Advection
In this section, it is assumed that, although the hydraulic flow though the orifice strongly affects the chemical flow, the chemical potential has negligible effect on the hydraulic pressure across the orifice. It is therefore possible to approximate the hydraulic dynamics by neglecting the chemical dynamics and, as discussed in this section, the chemical dynamics can be reformulated in terms of a modulated resistance, thus bringing transport of a chemical via advection within the same conceptual framework as transformation via a chemical reaction. The bond graph R:r i component represents the orifice hydraulic resistance, and the bond graph modulated transformer component TF:c i transduces energy between the chemical and hydraulic domains with modulus upstream concentration c i . The bond direction corresponds to energy flow from chemical to hydraulic and P i+1 = P i + ΔP i − Δ R P i and μ i+1 = μ i − Δμ i . Note that the total pressure drop P i − P i+1 might not be the same sign as the total chemical potential drop μ i − μ i+1 due to the pressure drop Δ R P i across the hydraulic R:r i component.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220492 If the hydraulic flow Q i is positive, the advective flow f i is given by equation (2.1); if hydraulic flow Q i is negative, the advective flow f i is dependent on the concentration c i+1 rather than c i . Moreover, the chemical potential μ i is given in [9] where K is a species-dependent constant [9]. Hence the advective flow f is given by 2) can be modelled as the constitutive relation of a bond graph modulated (by hydraulic flow Q i ) resistive R component; to reflect the special properties arising from equation (3.2) this is given a special name: the RA (advective R) component. This component is compared and contrasted with the chemical transformation Re component [9] in §4.
The energy flow, or power P RA dissipated by the advective R component RA is the difference between the product of the flow and effort variables on the two bonds and equation (2.4). Note that since the transport of chemical is driven by the hydraulic flow, P RA may be negative, in contrast to a typical Re component which always has nonnegative power dissipation.

Flow through a pipe
Although the flow of fluid though a pipe can be modelled by a PDE, it is convenient to use a lumped model instead; this has been examined in the context of blood flow by Safaei et al. [17]. Figure 3a shows a lumped model where the pipe is divided into a number of compartments, and figure 3b gives the corresponding bond graph. The hydraulic-chemical interaction is one-way, via the inter-compartmental hydraulic flows Q i ; thus the chemical part of figure 3b does not depend on the details of the hydraulic model; only the flows are required. Thus, for example, inertial terms could be added to the hydraulic model by appending bond graph I components to the 1 junctions [17,24]. Advection via a pipe introduces a time delay into the system dynamics, and the lumped approximation of a pipe introduces an approximate pure time delay. In particular, if the pipe has N compartments, the total pipe volume V is In the case of identical compartments In the case of constant flow Q, the molar content V i c i of the ith compartment has, using (3.2), a rate of change Using transfer function analysis and thus replacing d/dt c i by sc i where s is the Laplace variable, it follows that 2)) giving one-way interaction from the hydraulic domain to the chemical domain to represent advection of the chemical species by the fluid.

Coupled advection and transformation
As mentioned in the Introduction, there are two distinct meanings of 'molar flow': the transport of a substance from one place to another via advection and the transformation of one substance to another via a chemical reaction [10]. As discussed in §3, advection may be modelled using the RA bond graph component together with the capacitive Ce components and bonds carrying chemical potential μ and molar flow f. Moreover, transformation can be modelled by the Re bond graph component together with the capacitive Ce components and bonds carrying chemical potential μ and molar flow v [9]. Hence both concepts can be combined in a single bond graph; this combination is illustrated using a simple example. Further, it is shown how transport via the RA component and transformation via the Re component can be analysed within a common framework. For example, figure 5a-b corresponds to a substance A transformed to substance C advected though an orifice with hydraulic flow Q into a well-stirred compartment within which the substance C is transformed to substance B. The substance C occurs at two locations so the corresponding Ce components are distinguished by denoting them Ce:C 1 and Ce:C 2 ; the corresponding amounts of substance (C) in each location are x 1 and x 2 with corresponding concentrations c 1 and c 2 .
In the case of positive flow, and in the context of figure 5b, equation (3.2) becomes where f is the advective flow though RA:r and c 1 is the concentration corresponding to Ce:C 1 . Figure 5c is identical to figure 5b except that RA:r is replaced by Re:r. In this case, the transformation flow v though Re:r is [9, eqn. 2.2]: κ 0 is the reaction rate-constant, and m É C is the standard potential for the species. Equation (4.1) can be rewritten as where In other words, the RA component can be modelled as a Re component with flow-dependent rate-constant κ 0 and a one-way reaction. Equations (4.2) and (4.4) can be combined using the symbol λ where λ = 0 corresponds to equation (4.4) for advection and λ = 1 corresponds to equation (4.2) for transformation; equations (4.2) and (4.6) can be put in a common format as   Figure 6 illustrates the differences between the systems of (b,c).

Illustrative example
To illustrate the difference between transportation (λ = 0) and transformation (λ = 1), consider the special case of the system of figure 5 where the rate constants of the two reactions r 1 and r 2 are the same: κ 1 = κ 2 .
In the steady state, the time derivatives are zero and so, using equation (4.10) ð4:11Þ If the flow Q is large, κ is also large and equation (4.8) implies that c 1 % lc 2 : ð4:12Þ Hence large flow in the steady state implies that ð4:13Þ Thus when λ = 0 (advection, figure 5b)

ð4:14Þ
and when λ = 1 (transformation, figure 5c). ð4:15Þ The behaviour for a particular set of values is illustrated by the simulation results of figure 6. The steady-state values for large flow are explained by the analysis leading to equations (4.14) and (4.15).

Circulatory advection
The cardiovascular systems of humans are closed in the sense that the blood is recirculated round the body: it forms a circulatory advection system where the blood advects many substances, notably oxygen (O 2 ) and carbon dioxide (CO 2 ). There are, in fact, a number of circulatory systems including coronary circulation, pulmonary circulation, cerebral circulation and renal circulation. The transport of O 2 involves the binding and unbinding of O 2 to haemoglobin which is itself advected by the blood.
Using a simple single circulation model, this paper shows how, in principle, the advection models of §3.1 and the advection-transformation models of §4 can be used to build models of a circulatory advection system involving binding and unbinding. More complex systems could be built using the modular properties [6] of bond graphs. Figure 8 shows a simplified binding/unbinding cycle which has been split into two halves corresponding to two different locations. This section looks at ways in which these two halves can be connected using the advection models developed in §3. In each case, four advection connections are required to carry both bound C and unbound E enzymes to and from the two locations. Figure 9a shows orifice advection connections using four instances of the RA component of §3 whereas figure 9b uses four instances of the Pipe component of § 3.1.
If the forward and reverse orifices are identical, the two pairs of RA components in figure 9a, can be replaced by two Re components as shown in figure 9c; this is not true if the orifices are different or if pipe components are involved. Figure 10 gives some simulation results with unit parameter values except for volume V = 5 and (in the case of the pipe connection) number of lumps N = 5.

Pharmacokinetics
Pharmacokinetics is the study of drug uptake in living creatures; for example, the anaesthetic gas nitrous oxide (N 2 O) is administered during operations by adding it to inspired air. As it acts on the brain, the dynamics of the transfer of N 2 O from the lungs to the brain via blood circulation is of interest. In his classic papers, Mapleson [18,19] gave a fully parametrized compartmental model of the uptake of the anaesthetic gas N 2 O in humans outlined in figure 11. The  Figure 8. Split binding/unbinding cycle. At location 1, the small molecule A binds to enzyme E to form complex C; at location 2, the complex C unbinds into the small molecule A and enzyme E. C and E are transferred back and forth between the two locations by advection. For example, the small molecule A could represent oxygen and C and E could represent haemoglobin with and without bound oxygen. Location 1 could be the lungs and location 2 an oxygen-consuming organ.

Re:r1
Re  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220492 eight compartments are listed in figure 11; arterial blood flows from the lungs to the seven organ compartments and returns as venous blood. A simple bond graph model was given by Worship [25] and Gawthrop & Smith [26,ch. 9]. Rather than model the forking arteries and veins, following [19], the blood flow is approximated by a separate arterial flow to, and venous flow from, each compartment; the fraction δ p of blood passing into each compartment is listed in table 1. In this section, the bond graph model of [26] is extended to explicitly include the dynamics of the veins and arteries using the approach of §3.1. For each organ, the compartment Ce component has a constant K p given by As discussed in §3.1 and with reference to table 1, the pipe parameters are The modular bond graph model of figure 11 was simulated and the results appear in figure 12.
The results of [19, figure 2] for the three pools (kidney, grey matter and white matter) were digitized from an image of the printed paper using the online tool https:// automeris.io/WebPlotDigitizer/. As the digitized data have units of mmHg, they were normalized to be commensurate with the simulation and are shown as dots in figure 12a. The simulation and digitized results of figure 12a correspond closely; the discrepancies could be due to digitization errors, the more detailed pipe model used here, or the use of a different integration method from that used, but unspecified, by Mapleson [19].
The results shown in figure 12b indicate that N 2 O is stored in the fat; the slow release of N 2 O from the fat into the bloodstream can lead to unwanted post-operative anaesthesia [27]. In particular, at the conclusion of surgery, obese patients may wake more slowly than lean patients.

Conclusion
We have shown how advection can be incorporated into bond graph models of physiological systems, providing a physically consistent framework for linking advective flows with the fluxes associated with chemical reactions. This approach could be used to understand the function of a wide range of organs. Some examples include the transport of oxygen in the blood via haemoglobin; the transport of nutrients and drugs through the liver; and the transport of gases through the lungs. We believe that our approach will allow the incorporation of chemical transformation in these systems as well as transport [17,28,29]. Such models could find use in clinical contexts; for example predicting the functional consequences of surgically removing parts of the liver for cancer treatment [30] and understanding the effects of pulmonary obstruction on lung function [31,32].
Because bond graphs have already been used to model biochemical reactions [9], the models in this paper could be extended to incorporate tissue metabolism as well as advection. This has potential applications in studying the processing of nutrients by the liver and its consumption by tissue [33], as well as in understanding the metabolism royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220492 and clearance of drugs [34,35]. The approach could also be used to model complex biochemical kinetics, for example, the cooperative binding/unbinding of oxygen to haemoglobin [36]. A benefit of the approach in this paper is that the building blocks are modular. If a finer level of detail is required, one could replace the simple models of circulation in figure 11 with more realistic models of vasculature. Furthermore, several organs have a repeating hierarchical structure that can be exploited by a modular approach. For example, the lung is composed of branching pipes that conduct the flow of gas into millions of alveoli [31]. Similarly, the liver is composed of lobes which are in turn made up of several liver lobules [30]. Once bond graph models of fundamental building blocks have been developed, they can be re-used in constructing models of more complex anatomical structures [5]. Such an approach could prove valuable for developing the multi-scale models required for the Physiome Project.
To improve the realism of the models presented here, there are some issues that require further investigation.   Figure 11. Pharmacokinetics: (a) the schematic diagram shows how the lungs are connected to six compartments representing organs together with a shunt representing blood bypassing organs [18,19]. (b) A modular bond graph representation of (a). The pool module appears in (c) and seven instances appear.    [19], N 2 O is added to inspired air for the first 2 min; for the rest of the simulation pure air is inspired. (a) The concentration of N 2 O, normalized by the initial concentration in the inspired air, is plotted for the alveolar blood and the kidney, brain grey matter and white matter compartments. As can be seen, the N 2 O diffuses into and out of the blood via the alveoli and the concentration rises for the first 2 min and then falls as N 2 O is withdrawn. Having passed though the arteries, the arterial N 2 O then drives the individual compartments and is removed by the venous blood. The results of [19, figure 2] (digitized from an image of the printed paper and normalized to be commensurate with the simulation) are shown as dots. The effective time-delay of the arteries is t ¼ V A =Q % 13 s ¼ 0:22 min; this delay is reflected in the initial response of the three compartmental concentrations. (b) N 2 O accumulates in the fat and takes a long time to decay (beyond the range of this simulation).
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220492 1. Flow typically varies across the transporting medium and thus three spacial dimensions are required for a full description of advection to include, for example, shearinduced Taylor diffusion and the effect of the Womersley number. The bond graph approach can be used to model three-dimensional fluid flow [37]; however, as discussed by Safaei et al. [16,17], practical full-body haemodynamic models require both one-and three-dimensional flow models; moreover, one-dimensional models can provide boundary conditions for three-dimensional models [38]. This paper focuses on modelling advection in the onedimensional case using lumped (zero-dimensional) approximations, but future work could be directed to including advection within the three-dimensional case. The Port-Hamiltonian approach [39] provides one way forward. 2. We have assumed that the chemical domain has negligible effect on the hydraulic domain ( §3). While this is in line with prior approaches [40], the validity of this assumption needs further investigation. 3. This paper considers incompressible flow through rigid pipes; in the case of blood flow the pipes are not rigid, requiring extension to compartments with variable volumes; this will be the subject of future research. 4. Hydraulic flow reversal leads to the switching function of equation (3.2). This requires further analysis, possibly using switched bond graph methods [41].
5. Species may be carried in cells within the fluid stream-for example, haemoglobin is carried by red blood cells, which may require the modelling of separate compartments.
Data accessibility. The figures and tables in this paper were generated using the Jupyter notebooks and Python code available at https:// github.com/gawthrop/Advection22.