Energy-based Modelling of the Feedback Control of Biomolecular Systems with Cyclic Flow Modulation
EEnergy-based Feedback Control of BiomolecularSystems with Cyclic Flow Modulation
Peter J. Gawthrop *1,21
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, Universityof Melbourne University of Melbourne, Victoria 3010July 30, 2020
Abstract
Energy-based modelling brings engineering insight to the understanding of biomolecularsystems. It is shown how well-established control engineering concepts, such as loop-gain,arise from energy feedback loops and are therefore amenable to control engineering insight.The approach is illustrated using a class of metabolic cycles with activation and inhibitionleading the concept of Cyclic Flow Modulation.
The bond graph implementation of Network Thermodynamics was introduced some 50 years agoas an energy-based approach to modelling biomolecular systems [1]. “Graphical representationssimilar to engineering circuit diagrams can be constructed for thermodynamic systems. ... suchdiagrams do increase one’s intuition about system behaviour.”[2].The design of linear feedback circuits also has a long history and the correspondingly well-established theory of control systems [3] has been applied to biomolecular systems [4–6] andhas led to a number of control concepts such as feedback and integral action being used in thebiomolecular context [5].Classical linear control theory is based on transfer function models of dynamical systems. Incontrast, the energy-based approach of this paper uses the bond graph paradigm for modelling * Corresponding author. [email protected] a r X i v : . [ q - b i o . M N ] J u l iomolecular systems. There has been limited work on the bond graph approach to control [7–12]. For this reason, a novel method is introduced to allow the transfer function based approachof classical linear control to be utilised in the analysis of feedback systems modelled by bondgraphs and thus combine energy-based modelling with control systems analysis.As discussed by Gawthrop and Crampin [13], the bond graph approach gives the set of non-linear ordinary differential equations describing the biomolecular system being modelled. Lin-earisation of non-linear systems is a standard technique in control engineering. Linearisationin the context of bond graph models of biomolecular systems was introduced by Gawthrop andCrampin [13] and is used here.The role of metabolic cycles in the regulation of metabolic flux is well established [14–17].Such cycles are involved in a number of substrate conversions including those between fructose-6-phosphate and fructose-1,6-biphosphate, fructose-6-phosphate and fructose-2,6-biphosphate,triglyceride/fatty acid, glucose and glucose-6-phosphate, and glycogen and glucose 1-phosphate[14, 16, 17]. To illustrate the fusion of network thermodynamics and control theory, this paperwill focus on the first two inter-conversions involving fructose-6-phosphate (F P). Because of thecyclic nature of these two reactions, and the fact that flow is modulated, the term
Cyclic FlowModulation (CFM) is used to describe such reaction systems.The use of CFM requires energy and there is a trade-off between quality of control andenergy consumed [15]. It is therefore important to account for energy flows when modellingbiomolecular systems and this is done here using the fusion of the network thermodynamicsparadigm, as implemented using bond graphs, with control theory.Criteria for robust biochemical reaction networks have been established which ensure zerosteady-state error [18–20]; but these papers make no mention of energy and therefore entirelyignore thermodynamic constraints and the consequent performance-energy trade-off.Building complex systems is simplified using modularity; but it is essential to distinguishtwo different concepts of modularity: computational modularity where physical correctness isretained and behavioural modularity where module behaviour (such as ultra-sensitivity) is re-tained [13]. As well as providing computational modularity, bond graphs provide a natural for-mulation of behavioural modularity and reveal the sources of retroactivity [13].
Chemostats [13, 21] are used to create an open system from a closed system and also provide a convenientway of providing ports to connect bond graph modules.§ 2 introduces the bond graph based approach to the analysis of feedback control systemsusing an enzyme catalysed reaction with competitive inhibition as an illustrative example. § 3shows how cyclic flow modulation (CFM) can be used to build effective feedback controllerswith approximate integral action. § 4 concludes the paper and gives directions for future work.
Figure 1(a) depicts a conventional feedback control system in transfer-function form. The fourtransfer functions G con ( s ) , G sys ( s ) , G w ( s ) and G g ( s ) represent the controller, the system un- The pejorative term “futile cycle” is often used to describe such cycles; this will be avoided in this paper. − + + G con G sys yG d G w w ud (a) Transfer function [B] [P][S][Inh][Act] (b) Bond graph Figure 1: Feedback control. (a) A classical feedback loop block diagram representing the lin-earisation of a non-linear system. The blocks represent transfer functions which are connectedby signals. y is the controlled output, d a disturbance and w the setpoint, or desired value of y .(b) A bond graph feedback loop. CON and
SYS are bond graph modules with ports denotedby []. Ce : P , Ce : D and Ce : P0 are bond graph components representing species correspondingto product, disturbance and reference species respectively. The (cid:43) symbol indicates an ener-getic connection between two subsystems; the half-arrow indicates the direction correspondingto positive energy flow. 3er control, the setpoint and disturbance transfer functions respectively where s is the Laplacevariable. The four signals y , u , w and d represent the system output, system input, setpoint anddisturbance respectively.The closed-loop transfer function is: y = L ( s )1 + L ( s ) G w ( s ) w + 11 + L ( s ) G d ( s ) d (1)where L ( s ) = G con ( s ) G sys ( s ) (2) L ( s ) is referred to as the feedback loop gain . In the engineering context, G con ( s ) and G sys ( s ) would arise from separate physical entities; nevertheless, the loop gain L ( s ) (2) appearing inequation (1) only requires the product of G con ( s ) and G sys ( s ) . This is important for biomolecularsystems where there is no clear physical distinction between controller and system: it is thefeedback loop itself that is of fundamental importance.Typically, such control systems are analysed in the frequency domain by setting s = jω where j = √− and ω is frequency in rad / sec . At those frequencies where L ( jω ) is large,equation (1) can be approximated by y ≈ G w ( s ) w . In other words, a large loop gain L ( jω ) isdesirable insofar as the system output y is a close match to the desired value G w ( s ) w despitedisturbances represented by d . However, incorrect choice of the the loop gain L ( s ) can lead toinstability and L ( s ) is, moreover, subject to fundamental constraints [3].To summarise, there are two potentially conflicting issues in controller design: good distur-bance rejection and stability; these are both captured in the loop gain L ( s ) .Figure 1(a) implicitly assumes that the connection between subsystems, such as those rep-resented by G con ( s ) and G sys ( s ) is one-way as indicated by the arrows. However, the physicalcontroller needs to be designed to make sure this one-way interaction is correct; this requiresthe use of energy. It has been argued that this approach is misguided, even in the context ofengineering systems. This has lead to the concept of physical-model based control [7–12].In the context of biomolecular systems, the concept of retroactivity [6] has been introducedto explain why interaction is not one-way and thus design based on simplistic application of theapproach of Figure 1(a) often fails.There are two reasons why the bond graph approach is superior to the transfer function ap-proach of Figure 1(a) in the context of feedback control:1. It explicitly accounts for the two-way interaction found in physical systems in general andbiomolecular systems in particular.2. It explicitly accounts for energy flows and thus can directly expose performance/energyconsumption trade-offs.For this reason, the transfer function paradigm of Figure 1(a) is replaced by the bond graph basedparadigm of Figure 1(b).Figure 1(b) is based on the notation for modular bond graphs [22]. The two bond graph mod-ules are CON and
SYS ; CON represents the controller and has three ports: [Act] (activation),[Inh] (inhibition) and [Con] (control signal) and
SYS represents the system and has two ports:4S] (substrate) and [P] product. In the sequel, the controller module will be instantiated by threemodules in turn: an enzyme catalysed reaction with competitive activation and inhibition (§ 2.2),cyclic flow modulation (§ 3) and cyclic flow modulation with integral action (§ 3.1).The (cid:43) symbol indicates an energetic connection between two subsystems; the half-arrowindicates the direction corresponding to positive energy flow. In the biomolecular context, eachsuch bond is associated with two covariables: chemical potential µ J mol − and flow v mol / sec .The key point is that the product of µ and v is power p = µv W . Alternatively, it is possible toscale these co-variable by Faraday’s constant F C mol − to give φ = F µ V and f = F v A where J C − has been replaced by the more convenient unit V and C / sec has been replaced by the moreconvenient unit A [23]. The components Ce:P and
Ce:S represent the product and substratespecies respectively and the components
Ce:P0 and
Ce:D represent the reference species andproduct disturbance respectively; because
Ce:P0 and
Ce:D represent exogenous variables, theyare chemostats [21, 13].As shown in the sequel, the bond graph modelling approach can make use of the transferfunction approach to understand the dynamic properties of feedback systems of the form ofFigure 1(b). In particular, as shown in § 2.3, the fundamental control systems concept of loop-gain can be retrieved from the bond graph modelling paradigm. But first, linearisation must beconsidered.
Biomolecular systems are nonlinear and must be linearised before applying transfer functiontechniques. Linearisation of biomolecular systems in a biomolecular context, together with adiscussion on retroactivity, is given by [13]. In particular, the non-linear system equations are: ddt x = N f f = F ( x, x ch ) (3)In systems biology terms: the n x vector x represents the amount of each non-chemostattedspecies ( mol ), the n x × n f matrix N is the system stoichiometric matrix , the n f vector f rep-resents the flow in each reaction ( mol s − ).The n x ch vector x ch represents the amount of eachchemostatted species ( mol ). F ( x, x ch ) is a nonlinear function of both arguments. Because ofthermodynamic constraints, F has a particular structure dependent on the stoichiometric ma-trix N [24, 25] and is automatically generated from the bond graph representation. In standardcontrol system terms, x is the system state , f is the system output and x ch the system input .The corresponding linearised equations are: ddt ˜ x = N ˜ f (4) ˜ f = C ˜ x + D ˜ x ch (5)where ˜ x = x − ¯ x (6) ˜ x ch = x ch − ¯ x ch (7) ˜ f = f − ¯ f (8)5here the n f × n x matrix C and the n f × n x ch matrix D are given by the partial derivatives: C = ∂f∂x D = ∂f∂x ch (9)evaluated at the steady-state values ¯ x and ¯ f of state and flow respectively corresponding to theconstant chemostat state x ch = ¯ x ch : N ¯ f = 0 ¯ f = F (¯ x, ¯ x ch ) (10)Linearisation has two steps: finding the steady-state state ¯ x and flow ¯ f and then computing thelinearisation matrices C and D . The first is simply accomplished by numerically simulatingthe system until a steady-state is reached ( ddt x ≈ ). The second is achieved symbolicallywithin BondGraphTools ( https://pypi.org/project/BondGraphTools ) usingthe symbolic derivative functions of the sympy library ( ). ThePython Control Systems Library ( https://pypi.org/project/control/ ) is used toconvert the linearised system from state-space form to transfer function form, manipulate trans-fer functions and to generate time and frequency responses. The modified enzyme-catalysed reaction module of Figure 2(b) and the pathway module of Fig-ure 2(c) are embedded in the feedback loop of Figure 1(b) and used for the purposes of illustra-tion; the parameters are given in Figure 3.The non-linear system equations were derived from the modular bond graph of Figure 1(b)using
BondGraphTools and simulated to give the steady-state condition corresponding to theparameters of Figure 3. The linearised equations were then extracted and the transfer functionrelating the disturbance ˜ x D to the product ˜ x P generated. The corresponding closed-loop stepresponse appears in Figure 3.However, simulation does not provide an explanation of why the steady-state is the particularvalue shown nor why the dynamics are as shown. The explanation is provided by the analysis ofthe following section. As discussed above, the loop-gain L ( s ) is a key transfer function in the classical control systemsanalysis of Figure 1(a). This section indicates how the loop-gain L ( s ) can be derived from thebond graph of Figure 1(b).The closed-loop system of Figure 1(b) includes two chemostats Ce : P0 and Ce : D whichmake the corresponding states x P and x D independent variables; the product state x P remainsa dependent variable which evolves with time as in Figure 3. To create an open loop system,the component Ce : P representing the product is also made a chemostat thus making x P anindependent variable. 6 R e : r Ce:A Ce:BCe:Inh 1Ce:Act Ce:E00Ce:F Ce:G1 0Ce:C 010 R e : r R e : r Ce:E0 (a) Enzyme-catalysed reaction (ECR) [A] [B][Inh] ecr:ecr [Act] (b) Cooperativity Ce:S Ce:I1 Ce:I2 Ce:P R e : r R e : r R e : r (c) Pathway Figure 2: Enzyme-catalysed reaction (ECR). The Bond Graph notation is: (cid:43) energy connection; Ce species; Re reaction; common potential connection; common flow connection [26]. (a)Enzyme-catalysed reaction with competitive activation and inhibition. C : A , C : B , C : E , C : C , C : Act & C : Inh represent the substrate, product, enzyme, enzyme-substrate complex, activationand inhibition respectively. C : F & C : G , provide the driving energy. Example species appearin Figure 7. (b) A simple model of cooperativity is included by specifying that 4 activation and4 inhibition species interact with the enzyme; this is achieved in a modular fashion. (c) Thecontrolled system of the module Path:sys of Figure 4 is, in this example, a simple path of 3reactions represented by Re : r1 – R : r3 with intermediate species C : I1 and C : I2 .7 t x P ( n o r m a li s e d ) LinearisedNon-linear (1)Non-linear (0.1) g D = 0.23 Figure 3: Non-linear and linearised closed-loop step response. The asymptotic value g D is in-dicated by a dashed line. The amplitude of the disturbance step is given for the non-linear sim-ulations and the resultant response is divided by the step amplitude. The normalised nonlinearresponse is close to the linear case for an amplitude of 0.1 and differs slightly for an amplitude of1.0. For the purposes of illustration, all parameters are unity except K F = 10 , K G = 10 − , κ con and κ sys = 10 The steady-state values were x S = 12 . and x P = 10 . . +−+ G D ( s ) 1 sG P ( s ) G P ( s ) x P x D x P v P (a) Loop analysis +−+ + + G D ( s ) 1 sG act ( s ) G pas ( s ) G P ( s ) x P x D x P v P (b) Split-loop analysis Figure 4: Loop analysis. (a) The block diagram corresponding to opening the bond graph feed-back loop by setting the product Ce : P to be a chemostat. x P is the amount of product and v P the product flow. G P ( s ) , G D ( s ) and G pP ( s ) are transfer functions relating x P , x D and x P to f P . (b) The transfer function G P ( s ) is split into two terms: G act ( s ) and G pas ( s ) correspondingto active and passive feedback. 8he linearised flow ˜ f P into the chemostat Ce : P is given by the sum of three terms corre-sponding to the three chemostats Ce : P , Ce : P0 and Ce : D respectively: ˜ f P = − G P ( s )˜ x P + G P ( s )˜ x P + G D ( s )˜ x D (11)where − G P , G P and G D are the transfer functions relating ˜ f P to ˜ x P , ˜ x P and ˜ x D respectively.The minus sign associated with G P is to give compatibility with standard definitions of loop gainin a negative feedback context.To reclose the loop, Ce : P is restored to non-chemostatted dynamics using the transfer func-tion relating ˜ x P to ˜ f P : ˜ x P = 1 s ˜ f P (12)The block diagram corresponding to Equations (11) and (12) is shown in Figure 4(a). UsingEquation (2), the loop gain L ( s ) is given by: L ( s ) = G p ( s ) s (13)From the block diagram of Figure 4(a), or from Equations (11) and (12), the closed-loopsystem can be explicitly written as ˜ x P = 1 s + G P ( s ) [ G P ( s )˜ x P + G D ( s )˜ x D ] (14)The steady state value ¯˜ x P of ˜ x p is obtained by setting s = 0 to give ¯˜ x P = 1 G P (0) [ G P (0)¯˜ x P + G D (0)¯˜ x D ] (15)In particular, the steady-state disturbance gain g D is given by: g D = ¯˜ x P ¯˜ x D = G D (0) G P (0) (16)With the parameters given in Figure 3: G D (0) = 1 . (17) G P (0) = 4 . (18)hence g D = 0 . (19)This corresponds to Figure 3. 9 SYS 0 Ce:PCON Ce:S Re:rdCe:DCe:P0 Ce:Inh [B] [P][S][Inh][Act]
Figure 5: Split-loop. The feedback bond in Figure 1(b) is removed and replaced by the chemostat Ce : Inh . This splits the loop and allows active and passive feedback to be distinguished asdescribed in the text.
The previous section shows how the loop gain L ( s ) may be derived from the closed-loop systemin the bond graph form of Figure 1(b). This section expands this analysis by dividing the loopgain L ( s ) into two parts: and active part L act ( s ) and a passive part L pas ( s ) so that L ( s ) = L act ( s ) + L pas ( s ) (20)The active part arises mainly from the properties of the controller ( CON ); the passive part arisesmainly from the properties of the system (
SYS ) appearing in the closed-loop bond graph ofFigure 1(b).The split-loop procedure is based on removing the feedback bond linking the controlled prod-uct Ce : P to the inhibition port ([Inh]) of the controller. This is depicted in Figure 5 where thebond has been removed and the chemostat Ce : Inh has been added. To focus on the loop gain,the chemostats Ce : P0 and Ce : D are held at the steady state values ( ˜ x D = ˜ x P = 0 ) for therest of this section. The linearised flows ˜ f P P into the chemostat Ce : P and ˜ f II into the chemo-stat Ce : I are each given by the sum of two terms corresponding to the two variable chemostats Ce : Inh and Ce : P respectively: ˜ f P P = − G P I ( s )˜ x Inh − G P P ( s )˜ x P (21) ˜ f II = − G II ( s )˜ x Inh − G IP ( s )˜ x P (22)When the split-loop is reconnected ˜ x Inh = ˜ x P (23)and ˜ f P = ˜ f P P + ˜ f II (24) = − [ G P I ( s ) + G P P ( s ) + G II ( s ) + G IP ( s )] ˜ x P (25) G P I ( s ) is the transfer function from the inhibition port of the controller to the product and is thus10he active part of the control. Hence the previous equation is rewritten as: ˜ f P = − [ G act ( s ) + G pas ( s )] ˜ x P (26)where G act ( s ) = G P I ( s ) (27) G pas ( s ) = G P P ( s ) + G II ( s ) + G IP ( s ) (28)Once again, the minus signs associated with G act and G pass are to give compatibility withstandard definitions of loop gain in a negative feedback context.To allow comparison with Equation (11), the transfer functions appearing Equation (21) areevaluated with the same steady states as those of the closed-loop system and, in addition, recon-nection of the split loop implies ¯ x inh = ¯ x P ˜ x inh = ˜ x P (29)Comparing Equations (11 and Equation (21), it follows that: G P = G act + G pas (30)Further, defining L act ( s ) = G act ( s ) s L pas ( s ) = G pas ( s ) s (31)Equation (20) follows from Equation (30). Thus the block-diagram of Figure 4(a) can be ex-panded to give the block-diagram of Figure 4(b).The conventional approach to feedback control in the engineering context would regard L pas as an unwanted artefact to be eliminated by correct design; similarly, in the life-sciences context, L pas would be regarded as due to retroactivity and therefore undesirable [6]. A theme of thispaper is that both these attitudes are incorrect. In the engineering context, using such interactionsto improve control are well established as physical-model based control [7–12]. In the systemsbiology context, this paper will show that L pas has a stabilising influence on the control system.The closed-loop system is given in terms of G P ( s ) by Equation (14). Because of the de-composition (30), it is possible to see how the control system would, in principle, behave withonly passive or only active control. In particular, if ˜ x pass andn ˜ x act are the product concentrationdeviations in the two cases: ˜ x pass = 1 s + G pass ( s ) [ G P ( s )˜ x P + G D ( s )˜ x D ] (32) ˜ x act = 1 s + G act ( s ) [ G P ( s )˜ x P + G D ( s )˜ x D ] (33)Figure 6(a) shows G P ( jω ) , G pass ( jω ) and G act ( jω ) plotted on a Bode diagram [3]. The magni-tude of G act ( jω ) is large at low frequencies and small at high frequencies whereas the magnitudeof G pass ( jω ) is small at low frequencies and high at high frequencies. Hence the G P ( jω ) is closeto G pass ( jω ) at high frequencies yet retains the high gain at low frequencies due to G act ( jω ) .11 M a g n i t u d e Frequency (rad/sec)180135904504590 P h a s e ( d e g ) G P G act G pas (a) Components of G P ( s ) L ( j )3210123 I m a g L ( j ) x P = 1.0 L ( j ) L pas ( j ) L act ( j )Unit circle (b) Open-loop frequency response t x P x P = 1.0 Passive+Active (0.23)Passive (1.0)Active (0.29) (c) Closed-loop step response
Figure 6: ECR: Split-loop analysis. (a) Components of G P ( jω ) . G P I ≈ and is not shown.Of the two remaining components of G pas , G II is small at low frequencies and so G pas ≈ G P P at low frequencies; at mid and high frequencies, G pas provides phase advance. G act is large atlow frequencies and small at high frequencies. Thus G P ( jω ) ≈ G act ( jω ) + G P P ( jω ) at lowfrequencies and G P ( jω ) ≈ G pas ( jω ) at high frequency. Thus G act ( jω ) provides high gain (andthus low steady-state error) at low frequencies and G pass ( jω ) provides stabilising phase advanceat mid frequencies. (b) Open-loop frequency response. The loop gain L ( jω ) = G P ( jω ) s and itstwo components L act ( jω ) and L pas ( jω ) are plotted on a Nyquist diagram. L act ( jω ) passes closeto the − point and thus, without the term L pas ( jω ) , L ( jω ) would correspond to closed-loopsystem close to instability. However, at the relevant frequencies, L ( jω ) ≈ L pass ( jω ) and is wellaway from the − point and thus corresponds to a stable closed-loop system. (c) Closed-loopdisturbance step response. The hypothetical closed-loop responses corresponding to L pass ( jω ) and L act ( jω ) are well damped with large steady-state error and oscillatory with small steady-stateerror respectively. The actual response corresponding to L ( jω ) combines the best of both.12n classical control theory, the frequency response of the loop gain reveals dynamical prop-erties – including stability – of the closed-loop system. In particular, s is replaced by jω where j = √− and ω is frequency in rad s − . One such frequency-based approach is based on theNyquist diagram [3] where the imaginary part of L ( jω ) is plotted against the real part of L ( jω ) for a range of frequencies. The phase ∠ L ( jω ) when the modulus | L ( jω ) | = 1 is of interest,hence the unit circle is plotted on the Nyquist diagram of Figure 6(b). There are three frequencyresponses plotted: L pass ( jω ) shows that, as the frequency response is well away from the − point, the time response is well-damped; L act ( jω ) shows that, as the frequency response passesclose to the − point, the time response is oscillatory; L ( jω ) is similar to L pass ( jω ) near the unitcircle and therefore also has a well-damped response.The corresponding unit step responses appear in Figure 6(c) along with the step response of ˜ x P corresponding to Equation (14). The disturbance response of the passive-only system is well-behaved but the steady-state value is large; in contrast, the disturbance response of the active-onlysystem is oscillatory but the steady-state value is small. The overall controller combines the bestof both responses: it is well behaved with a small steady-state value. The numerical steady-statevalues for the overall controller are given in Equation (19); in a similar fashion: g pass = 1 G pass (0) = 1 / . (34) g act = 1 G act (0) = 1 / .
44 = 0 . (35)(36)Thus the small steady-state value is largely due to the active part of the control. “ The parallel existence of two irreversible reactions is of the greatest importance in metabolicregulation: it means that the direction of flux between two metabolites is determined by differ-ential regulation of the activities of the two enzymes ” [16]. A bond graph interpretation of thismechanism appears in Figure 7(a) and this will be used as the basis replacing the
CON compo-nent in the bond graph feedback loop of Figure 1(b) by a more sophisticated control actuator.The use of such cyclic flow modulators is motivated by the pair of key metabolic reactionsdiscussed by Cornish-Bowden [16]:F P + ATP
PFK F P + ADP (37)F P + H O PBP F P + Pi (38)This pair of reactions can be related to the CFM bond graph of Figure 7(a) (with reference toFigure 2) as follows. The enzyme corresponding to
ECR : Fwd is PFK (phosphofructokinase)and the enzyme corresponding to
ECR : Rev is FBP (fructose biphosphatase). The substratecorresponding to Ce : A is F P (fructose-6-phosphate) and the product corresponding to Ce : A [B][A][Inh][Act][A][B] [Act][Inh] ECR:Rev (a) Cyclic flow modulation (CFM) C e : A c t
01 100Ce:A Ce:B11 0 C e : I n t C e : I nh [A] CFM:I [Inh][B][Act][A] [B][Inh]
CFM:P [Act] (b) Integral action (CFMI)
Figure 7: Cyclic flow modulation. (a) The two components
ECR : Fwd and
ECR : Rev are in-stances of the modulated ECR of Figure 2(b). The Ce : A component represents both the substrateof ECR : Fwd and the product of
ECR : Rev ; the Ce : B component represents both the substrateof ECR : Rev and the product of
ECR : Fwd . Component Ce : Act both activates
ECR : Fwd and inhibits
ECR : Rev ; component Ce : Inh both inhibits
ECR : Fwd and activates
ECR : Rev .(b)
CFM : P gives proportional (P) action whereas CFM : I gives integral (I) action by driving thespecies Ce : Int which activates
CFM : P . To exemplify strong activation, three activation bondsare used. 14s F P (fructose-1,6-biphosphate). Within
ECR : Fwd , Ce : F corresponds to ATP (Adenosinetriphosphate) and Ce : G corresponds to ADP (Adenosine diphosphate); within ECR : Rev , Ce : F corresponds to H O and Ce : G corresponds to Pi (inorganic phosphate).Species which activate PFK and inhibit FBP include AMP (Adenosine monophosphate) andF P (fructose-2,6-phosphate); species which inhibit PFK and activate FBP include ATP and Cit(citrate).This section examines the effect of replacing the ECR based control of the feedback loop ofFigure 1(b) by a CFM based controller. The ECR module of Figure 2 has four visible chemostats Ce : A , Ce : B , Ce : Act , and Ce : Inh the latter three of which are used as ports ([B], [Act], [Inh])in the feedback loop of Figure 1(b); the same chemostats ( Ce : A , Ce : B , Ce : Act , and Ce : Inh )are visible in the CFM module of Figure 7(a) and can be used as ports in the same way. Thus theCFM module can directly replace the ECR module in the feedback loop of Figure 1(b) (whichwas analysed in Figure 6); this CFM-based feedback loop is analysed in Figure 8.The linearised response of the ECR and CFM are similar for the parameters chosen. However,there is a significant difference: the CFM controller is bidirectional, the ECR is not. In bothcases, the constant low-frequency gain corresponds to the proportional (P) controller of classicalcontrol. In contrast, the next section shows that two CFMs can be combined to give the classicalproportional+integral by endowing the controller with integral action . Integral action is an important concept in classical control theory [3] and endows a control systemwith zero steady-state error. In section of their paper Cloutier and Wellstead [5] discuss the role of F P (fructose-2,6-biphosphate), a strong activator of PFK (phosphofructokinase). In particular, F P interconvertswith F P (fructose-6-phosphate) via the reaction cycle:F P + ATP
PFK2 F P + ADP (39)F P + H O F26BP F P + Pi (40)catalysed by the enzymes PFK2 (phosphofructokinase-2) and F BP (fructose-2,6-biphosphatase).The species which simultaneously activate PFK2 and inhibit F26BP include AMP and F P.Hence this pair of reactions is a further example of Cyclic Flow Modulation (CFM).Moreover, the PFK CFM and the PFK2 CFM strongly interact: the PFK CFM is positivelymodulated by the product of the PFK2 CFM: F P and both are positively modulated by AMP.Figure 7(b) gives the bond graph abstraction of the two interacting cycles.
CFM : P corre-sponds to the PFK-based CFM giving proportional (P) action whereas CFM : I corresponds tothe PFK2-based CFM and, as will be seen, gives integral (I) action. Within each CFM, the in-terpretation of the species is the same except that the product B of CFM : I corresponds to F Prather than F P. The component Ce : Int corresponds to the product F P which then activates
CFM : P . For illustration, and to emphasise the strong activation effect, three bonds represent theactivation effect. 15 M a g n i t u d e Frequency (rad/sec)1359045045 P h a s e ( d e g ) G P G act G pas (a) Components of G P ( s ) L ( j )3210123 I m a g L ( j ) x P = 1.0 L ( j ) L pas ( j ) L act ( j )Unit circle (b) Open-loop frequency response t x P x P = 1.0 Passive+Active (0.24)Passive (0.61)Active (0.39) (c) Closed-loop step response
Figure 8: CFM: Split-loop analysis. Detailed comments and parameters are given in Figure 6;the low frequency gain is higher leading to a lower steady-state error. Once again, the passiveterm G pas stabilises the high-gain active control term G act .16 M a g n i t u d e Frequency (rad/sec)9045045 P h a s e ( d e g ) G P G act G pas (a) Components of G P ( s ) L ( j )3210123 I m a g L ( j ) x P = 1.0 L ( j ) L pas ( j ) L act ( j )Unit circle (b) Open-loop frequency response t x P x P = 1.0 Passive+Active (0.029)Passive (0.61)Active (0.03) (c) Closed-loop step response
Figure 9: CFMI: Split-loop analysis. Detailed comments and parameters are given in Figure 6.Compared to the CFM controller response of Figure 8, the low-frequency gain of the active term G act rises as frequency decreases; this is the behaviour expected of an integrator. However, as theintegrator is not perfect, the gain is not infinite at ω = 0 . This approximate integral effect givesa lower steady-state error than the CFM controller whilst the passive term G pas continues to actto give a damped response. As the phase of L act ( jω c ) is below − ◦ at the critical frequency ω c at which magnitude | L act ( jω c ) | = 1 , the closed-loop system corresponding to the active partof the controller is unstable. 17
10 20 30 40 50 t x P ECRCFMCFMI g D = 0.23 g D = 0.24 g D = 0.029 Figure 10: Controller comparison. The disturbance response of the three controllers: ECR, CFMand CFMI is shown: the steady state disturbance gains g D (16) are shown. Thus the CFMIcontroller is the best of the three in reducing the steady-state error.This section examines replacing the CFM based control within the feedback loop of Figure1(b) by a CFMI based controller. As in the case of CFM, the same chemostats as ECR arevisible in the CFMI module of Figure 7(b) and can be used as ports in the same way. Thus theCFMI module can directly replace the ECR module in the feedback loop of Figure 1(b) analysedin Figure 6; this CFMI-based feedback loop is analysed in Figure 9. Compared to the CFMcontroller response of Figure 8, the low-frequency gain of the active term G act rises as frequencydecreases; this is the behaviour expected of an integrator. However, as the integrator is notperfect, the gain is not infinite at ω = 0 ; but this approximate integral effect gives a significantlylower steady-state error than the CFM controller whilst the passive term G pas continues to act togive a damped response. The disturbance response of the three controllers in shown in Figure 10;the CFMI controller has a substantially smaller steady-state disturbance error than the other two. In the examples so far, the activation chemostat of Figure 1(b) is defined by a unit state x P = 1 .By analogy with the classical feedback loop of Figure 1(a), it would be expected that x P wouldplay a similar role to w . Figure 11(a) indicates that this is approximately true for the CFMIcontrol: ¯ x P ≈ x P . Furthermore, varying x P changes the steady-state product flow. In thiscase, as the disturbance reaction gain is κ rd = 1 the product flow ¯ f P = x P − x d = x P − . Oneof the benefits noted for CFM control at the beginning of § 3 is that bidirectional product flow ispossible: Figure 11(b) illustrates this for the CFM and CFMI controller; it is not possible for the18 x P x P CFMCFMI (a) ¯ x P x P f CFMCFMI (b) ¯ f P Figure 11: Steady-state. (a) The steady-state product state ¯ x P ≈ x P for the CFM and CFMIcontrollers. (b) The steady-state product flow ¯ f P is bidirectional in the case of the two CFM-based controllers; this is not possible for the ECR controller.ECR controller. 19 Conclusion
Network thermodynamic modelling via bond graphs has been amalgamated with classical controltheory. The dual roles of active and passive feedback have been analysed: active feedback givesgood steady state performance whereas passive feedback provides stabilisation.Cyclic flow modulation (CFM) has been motivated by the phosphofructokinase-fructosebiphosphatase reaction metabolic cycle and shown to have a modular bond graph representa-tion. CFM can be used to build the proportional (P) and proportional+integral (PI) controllers ofclassical control theory, as well as allowing bidirectional flow modulation.Future work will include building an energy-based model of metabolism with AMP feedbackand mitochondrial transduction using existing energy-based models [23, 27].An important potential result of combining control theory with energy-based modelling is toidentify performance/energy trade-offs. This is important to both evolutionary theory [28] andsynthetic biology [29].
Peter Gawthrop would like to thank the Melbourne School of Engineering for its support via aProfessorial Fellowship, and Edmund Crampin and Michael Pan for help, advice and encourage-ment.
References [1] George Oster, Alan Perelson, and Aharon Katchalsky. Network thermodynamics.
Nature ,234:393–399, December 1971. doi:10.1038/234393a0.[2] A.S. Perelson. Network thermodynamics. an overview.
Biophysical Journal , 15(7):667 –685, 1975. doi:10.1016/S0006-3495(75)85847-4.[3] Karl Johan Astr¨om and Richard M Murray.
Feedback systems: an introduction for scientistsand engineers . Princeton University Press, 2008. ISBN 978-0-691-13576-2.[4] Michael A. Savageau.
Biochemical Systems Analysis. A Study of Function and Design inMolecular Biology . Addison-Wesley, Reading, Mass., 40th anniversary issue edition, 2009.[5] Mathieu Cloutier and Peter Wellstead. The control systems structures of en-ergy metabolism.
Journal of The Royal Society Interface , 7(45):651–665, 2010.doi:10.1098/rsif.2009.0371.[6] Domitilla Del Vecchio and Richard M Murray.
Biomolecular Feedback Systems . PrincetonUniversity Press, 2014. ISBN 0691161534.[7] Dean Karnopp. Bond graphs in control: Physical state variables and observers.
Journal ofthe Franklin Institute , 308(3):219 – 234, 1979. doi:10.1016/0016-0032(79)90114-5.208] A. Sharon, N. Hogan, and D. E. Hardt. Controller design in the physical domain.
Journalof the Franklin Institute , 328(5):697–721, 1991.[9] P. J. Gawthrop. Physical model-based control: A bond graph approach.
Journal of theFranklin Institute , 332B(3):285–305, 1995. doi:10.1016/0016-0032(95)00044-5.[10] P.J. Gawthrop, D.J. Wagg, and S.A. Neild. Bond graph based control and substruc-turing.
Simulation Modelling Practice and Theory , 17(1):211–227, January 2009.doi:10.1016/j.simpat.2007.10.005.[11] P.J. Gawthrop, B. Bhikkaji, and S.O.R. Moheimani. Physical-model-based control of apiezoelectric tube for nano-scale positioning applications.
Mechatronics , 20(1):74 – 84,February 2010. doi:10.1016/j.mechatronics.2009.09.006.[12] Peter Gawthrop, S.A. Neild, and D.J. Wagg. Dynamically dual vibration absorbers: a bondgraph approach to vibration control.
Systems Science and Control Engineering , 3(1):113–128, 2015. doi:10.1080/21642583.2014.991458.[13] P. J. Gawthrop and E. J. Crampin. Modular bond-graph modelling and analysis of biomolec-ular systems.
IET Systems Biology , 10(5):187–201, October 2016. doi:10.1049/iet-syb.2015.0083.[14] E.A. Newsholme, R.A.J. Challiss, and B. Crabtree. Substrate cycles: their role in improvingsensitivity in metabolic control.
Trends in Biochemical Sciences , 9(6):277 – 280, 1984.doi:10.1016/0968-0004(84)90165-8.[15] H. Qian and D. A. Beard. Metabolic futile cycles and their functions: a systems analysisof energy and control.
IEE Proceedings - Systems Biology , 153(4):192–200, July 2006.doi:10.1049/ip-syb:20050086.[16] Athel Cornish-Bowden.
Fundamentals of enzyme kinetics . Wiley-Blackwell, London, 4thedition, 2013. ISBN 978-3-527-33074-4.[17] Reginald H. Garrett and Charles M. Grisham.
Biochemistry . Cengage Learning, Boston,MA, 6th edition, 2017.[18] Guy Shinar and Martin Feinberg. Design principles for robust biochemical reaction net-works: What works, what cannot work, and what might almost work.
Mathematical Bio-sciences , 231(1):39 – 48, 2011. doi:10.1016/j.mbs.2011.02.012. Special issue on biologicaldesign principles.[19] Stephanie K. Aoki, Gabriele Lillacci, Ankit Gupta, Armin Baumschlager, David Schwe-ingruber, and Mustafa Khammash. A universal biomolecular integral feedback controllerfor robust perfect adaptation.
Nature , 570(7762):533–537, 2019. doi:10.1038/s41586-019-1321-1. 2120] Jinsu Kim and German Enciso. Absolutely robust controllers for chemical reaction net-works.
Journal of The Royal Society Interface , 17(166), 2020. doi:10.1098/rsif.2020.0031.[21] Matteo Polettini and Massimiliano Esposito. Irreversible thermodynamics of open chemicalnetworks. I. Emergent cycles and broken conservation laws.
The Journal of ChemicalPhysics , 141(2):024117, 2014. doi:10.1063/1.4886396.[22] Peter J. Gawthrop, Joseph Cursons, and Edmund J. Crampin. Hierarchical bond graphmodelling of biochemical networks.
Proceedings of the Royal Society A: Mathematical,Physical and Engineering Sciences , 471(2184):1–23, 2015. doi:10.1098/rspa.2015.0642.[23] P. J. Gawthrop. Bond graph modeling of chemiosmotic biomolecular energytransduction.
IEEE Transactions on NanoBioscience , 16(3):177–188, April 2017.doi:10.1109/TNB.2017.2674683.[24] A. van der Schaft, S. Rao, and B. Jayawardhana. On the mathematical structure of balancedchemical reaction networks governed by mass action kinetics.
SIAM Journal on AppliedMathematics , 73(2):953–973, 2013. doi:10.1137/11085431X.[25] P. Gawthrop and E. J. Crampin. Bond graph representation of chemical reactionnetworks.
IEEE Transactions on NanoBioscience , 17(4):449–455, October 2018.doi:10.1109/TNB.2018.2876391.[26] Peter J. Gawthrop and Edmund J. Crampin. Energy-based analysis of biochemical cyclesusing bond graphs.
Proceedings of the Royal Society A: Mathematical, Physical and Engi-neering Science , 470(2171):1–25, 2014. doi:10.1098/rspa.2014.0459.[27] Peter J. Gawthrop, Peter Cudmore, and Edmund J. Crampin. Physically-plausiblemodelling of biomolecular systems: A simplified, energy-based model of the mito-chondrial electron transport chain.
Journal of Theoretical Biology , 493:110223, 2020.doi:10.1016/j.jtbi.2020.110223.[28] Nick Lane. How energy flow shapes cell evolution.
Current Biology , 30(10):R471 – R476,2020. doi:10.1016/j.cub.2020.03.055.[29] Hadrien Delattre, Jing Chen, Matthew J. Wade, and Orkun S. Soyer. Thermodynamicmodelling of synthetic communities predicts minimum free energy requirements for sulfatereduction and methanogenesis.