J. theor. Biof. (1979) 78, 4999517

Bond Graph Simulation of Skeletal Muscle Glucose Metabolism S. A. HUNTER, R. A. PEURA, T. C. AND R. J. HARVEY

CRUSBERC

Biomedical Engineering Program and Department qf’ L@ Sciences, Worcester Polytechnic Institute, Worcester, Muss. 01609, U.S.A. (Received 26 June 1976, and in revised, form 4 January 1979) A mathematicalmodel of skeletalmuscleglucosemetabolismis presented. This modelresultedfrom the application of thermodynamicscombinedwith the dynamics of an energy storage capability. The physiological system consistsof the insulin modulated,carrier-mediatedtransport of glucoseinto skeletal muscle, the biochemical reactions of the EmbdenMeyerhof pathway, the cyclic 3’,5’ adenosinemonophosphate-dependent protein kinasecontrolled synthesis/degradation of glycogenstores,and the diffusion of lactate from muscle.The metabolicsystemis definedand synthesizedby the construction of the energy flow bond graph. The bond graph model wasevaluatedby simulatingthe systemresponse over a 2 min period to step increasesin extracellular epinephrine concentration. The simulatedresponseof the metabolitesand modulated enzymes corresponds, qualitatively and quantitatively, with in rive measurements publishedin the literature.

1. Introduction

General systems theory as described initially by von Bertalanffy (1962) emphasizes the need to explain complex phenomena, such as an entire organism, through their organization and not through single entities of which they are composed. Until recently, most systems models have employed electrical network schematic notation (Incropera, 1975 ; Lake, 1967; Lommen et al., 1971) which include details of their chemical thermodynamics, interacting diffusive and convective flows of dissolved components, and discrete system elements. In particular, bond graphs (Karnopp & Rosenberg, 1975; Gebben, 1977) provide a system modelling technique which presents a topological graph of power and information interactions between components of a system. The bond graph permits direct derivation of differential equations describing system dynamics, through an analog model of a system. Auslander et al. (1972) and Oster et al. (1973) 499 0022.-5193/79/120499+

19 $02.00/O

IQ 1977 Academic

Press Inc. (London)

Ltd.

MGATP MGADP MG

-

“IDH+H+

GLllCOSE tnnsport

KREBS CYCLE

triase 1 phosphate isomerass

H+

N”+A;V-

PYf’J’qPYRUVATE

2,3-ryttme P-PtlOSvYCERATE

IJ-DIPHOSPHOGLYCERATE 3-phogphdglywde

LACTATE dif* EXTRACELLULAR

H20v~~~V~T~ MGADP+H+ MGATP

MGATP

H+

GLYCERA’ -,.iYDE-3 =’ .PM. .osPHATE m-3-

MGATP MGADP+

NAD-

MGADP+H+

NADH+

FRUCTOSE-&‘HDSPE -v-r FRUCTOSE-l,6 DIPHOSPHATE

ADENOSINE MONOPHOSPHATE odenylaje kinase ADENOSiNE DIPHOSPHATE MBADP

eNAD GLYCEROL-~-PHOSPHATE glycerol phosphate __ shuttle

“‘“““““YYACETONE H ;PHATE

GLYCOGEN metabolism

MTRACELLUAR insuh rnodubied

LACTATE

PI

(a)

BOND GRAPH

SIMULATION

OF METABOLISM

501

lb’

EPINEPHRINE !i!z%E TRIPHOSPHATE J// pyg!$g$TE

A;;-+

I ,o~$~~~$ATE-y

dE2i

y

*ZN&ZTE Ii+

k0

- 3’, 5’- ADENOSHE

R-CYCLIC

MOWPHOSPHATE INACTIVE PROTEIN KINASE

cl

/

ACTIVE ,PROTEIN / KINASE

‘\ ‘\

active

/I

r

/

’ protein kinase I ACTIVE PHOSPHORYLASE GLYCOGEN UDP

1

INACTIVE PHOSPHORYLASE KINASE \ \ I. \

PI

PYROPHOSPHATE UTP+H+

pnospnary!se

(P’asa

I

\

-II,

UDP

wcogen

i phosphatase

t,

\ , ,I

PI \

1 wrophosphwlase

GLucosE-

+I-I~SPHATE phosphogl~omutase GLUCOSE-6-PHOSPHATE

FIG. 1. Biochemical pathways, reactions and control points of mammalian skeletal muscle employed in deriving the physiological model. (a) Glycolysis and ancillary reactions. (b) Glycogen synthesis and degradation.

502

S. A.

HUNTER

ET

AL.

applied bond graph models to study network thermodynamics and biophysical systems. Allen (1978) has adopted bond graph formulation to describe certain biological and biophysical processes active in a portion of a growing plant. The various enzymatic reactions making up the EmbdenMeyerhof pathways of anaerobic glycolysis have been extensively studied, yet this sequence of biochemical events has never received the attention of a thorough systems analysis. Consideration of the several biochemical reactions of glycogen metabolism and the cyclic-3’,5’-adenosine monophosphate-dependent protein kinase control system magnifies the overall systems problem. The Embden -Meyerhof pathway together with glycogen metabolism in mammalian muscle seemed likely candidates for which to create a mathematical model, employing modern principles of systems theory. A recent report (Rapoport et al., 1976) described a simple mathematical model for glycolysis in erythrocytes using ordinary differential equations. The complexity of that model, in which alternate competing pathways and intracellular compartmentalization are absent, illustrates the potential magnitude of the problem to a systems approach employing more complicated cell types. This report describes the principles of mathematics and engineering employed while developing this predictive model for both energy and material flow in mammalian sarcoplasm, during conversion of glycogen to pyruvate. The specific reactions and transport phenomena comprising the physiological model are those of glycogen metabolism initiated by adenyl cyclase and proceeding to the formation of glucose-l-phosphate, and the conversion of glucose to lactate during glycolysis. Specific reactions and associated control points are illustrated in Fig. 1. Many assumptions concerning the behavior of the physiological system have, of course, been made. The physiological system was assumed to be buffered against changes in the concentration of hydrogen ion, inorganic phosphate, pyrophosphate, magnesium ion, uridine triphosphate and uridine diphosphate. Also, no distinction was made between the various configurations of the metabolites. All metabolites were, in addition, assumed to be homogeneously distributed throughout the cytoplasm. The metabolism of pyruvate via the reactions of the citric acid cycle and the transport of dihydroxyacetone phosphate and glycerol-3-phosphate between the cytosol and the mitochondrion via the glycerol phosphate shuttle were assumed to proceed at a constant velocity. The metabolism of glucose-6phosphate by reactions in the hexose monophosphate shunt was assumed to be negligible in comparison to other metabolic routes in skeletal muscle. The processes of lipolysis and proteolysis were assumed to remain unchanged and have no

BOND

GRAPH

SIMULATION

OF

METABOLISM

503

effect on changes in glucose metabolism. The initial state as defined by the metabolite concentration profile and reaction velocity profile obtained from several different animals species was assumed to correspond to that of an average skeletal muscle. The mathematical model of skeletal muscle glucose metabolism presented in this report is based upon a thermodynamic analysis of the energy transforming biochemical reactions coupled with the dynamics associated with an energy storage capability. The analysis resulted in a model which is analogous to a non-linear resistance-parameter modulated capacitance network. The resistive or energy dissipating elements correspond to the biochemical reactions and membrane transport phenomena, and the capacitive or energy storage elements result from the potential energy associated with material storage. The model is presented in the form of a bond graph, a diagram in which the flow of energy through a system is displayed. 2. Modelling (A)

MATHEMATICAL

Strategy FORMULATION

The thermodynamic analysis of the system consisting of three open compartments [muscle (M), membrane (D) and extracellular fluid (B)] illustrated in Fig. 2 resulted in an energetic model of skeletal muscle glucose metabolism. The investigation involved a study of each compartment operating as an independent system and a study of the three compartment system operating as an integral unit. From this analysis, the system was partitioned into dissipative elements and storage elements with the appropriate energy flow or power term interconnecting these elements. The skeletal muscle compartment (M) at a temperature T, and of electrostatic potential 1+5,contains several sequential reactions designated by the indices p (p = 1,. . ., x) and two interactions at the compartment boundary, d,Q,/dt and material transport dIniz/dt, where ni represents the number of moles of species i. The application of the principle of conservation of mass relates the rate of change of the moles of species i to the rate of material transport and to the chemical reaction velocity d&,/dt through the equation

(1) where viP is the stoichiometric coefficient of species i in the p reaction. The conservation of energy principle yields an expression balancing the rate of

504

S. A.

HUNTER

ET

AL.

I

Environment Extracellular fluid

Membrane

Muscle

i FIG. 2. Energy flow graph of mammalian skeletal muscle glucose metabolism describing the various components of the physiological model. Specific terms and notations are defined in the text.

change of the energy of the system (dE/dt), with the rate of the energetic transfers occurring at the compartment boundary. The modes of energy flow at this boundary include heat (d,Q,/dt), entropy flow coupled with material transport

and electrical power (Y,,Z) using hi to symbolize the specific enthalpy of species i and I to symbolize the time derivative of electrical charge. Conservation of energy results in the relationship hi d,niz +C~+w. I

(2)

BOND

GRAPH

SIMULATION

OF

METABOLISM

505

The thermodynamic analysis is based upon the postulate that the Gibbs equation, originally derived for the equilibrium state, gives a valid description of a system at the non-equilibrium state (Prigogine, 1967). For the non-equilibrium situation of the skeletal muscle compartment, the Gibbs equation is

with entropy designated by S and the chemical potential of species i represented by pi. Defining the chemical affinity of a reaction by A = -c Vipi, the time derivative of charge

where zi is the electrovalency of the ionic species i and .F the Faraday constant, and noting that charge is conserved in a chemical reaction, the combination of the conservation of mass equation (1) and the conservation of energy equation (2) with the Gibbs equation specifies the rate of change of entropy as _ 1

4Qz

1 (Pi-hi)

T,dt

d,ni,

i-,--z

Thermodynamic theory also postulates that the entropy of a system (a) is composed of an entropy flow interaction with the system’s environment d,S/dt and an entropy production due to irreversible processes contained with the system d,S/dt. The sum of the entropy flow and entropy production terms is equal to the rate of change of entropy of the system (4)

Considering the muscle compartment as a part of the overall system, the entropy flow to the environment is zero and the entropy production is

4s _ 1 c-4 -- T, dt

M

d,Q,

dt

c

(Pi-hi)

---dtT

dIni*

S. A. HUNTER

506

ET

AL

A similar analysis is used for both the membrane (D) and the extracellular fluid (B) compartments. The entropy of the entire system (a) is equal to the sum of the entropy of each compartment

The substitution

of equations for the appropriate zi.H’, --dt+’

i

(~i+~i.~Y),-(~i+zi,~Y), T, (~i+zi.~Y)B-(1Li+Zi,~y)o T7

Td

term in equation (6) yields dEni x’

(7)

d,niz dt d,nil dt

The entropy production (d,S/dt), defines the dissipative functions in the three-compartment model due to gradients existing in the chemical potential, electrical potential and chemical affinity. Dclfining the forward

ajfiniry

with 4 the stoichiometric

coefficient of reactant i and the reverse affinity

with v; the stoichiometric coefficient of the product i, the entropy production due to the chemical reaction may be rewritten as (9)

where &,-A’, is the affinity gradient. The graphical programming of equations (4), (8) and (9) in Fig. 2 resulted in an energy flow graph showing the energetic interactions between and

BOND

GRAPH

SIMULATION

OF

METABOLISM

507

within the three compartments associated with dissipative structures. The entropy production equation (8) identified the dissipative elements. The entropy conservation equation (4) required a summing junction, a symbol used to indicate that the output is the algebraic sum of the inputs. The arrows represent energy flow or power bonds. Bonds with a half arrowhead indicate the assumed direction of flow, but are reversible. Bonds with a full arrowhead represent unidirectional or irreversible energy flow dictated by the fact that the entropy production of a non-equilibrium process must be greater than zero. The effort or potential variable (intensive property) was written above or to the left of the bond. The flow or through variable (extensive property) was written under or to the right of the bond. The product of the conjugate effort-flow pair is the power transmitted. In this particular investigation the extracellular fluid and membrane compartments were assumed to remain in a steady-state, while the muscle compartment exhibited dynamics associated with an energy storage capability. The change in energy of this compartment was due to the change in internal energy U where U = TS-PV

+ C nisi + C =i.~Yni.

Using the assumption that the Gibbs equation is valid in the non-equilibrium state, the time derivative of the internal energy for the muscle compartment is

assuming a constant volume process. This equation indicates the modes of reversible energy storage and the corresponding effort-integrated flow pair. Energy storage which is a function of the effort-integrated flow pair is capacitive as opposed to inductive storage (Oster, Perelson & Katchalsky, 1973).

The storage flow arising from the change in mass, the chemical reaction flow and the material flow across the compartment boundary, are coupled through the conservation of mass equation

!-) dni

(12)

dt M Since the storage elements are reversible, conservation of energy applied to these elements results in the relationship (13)

508

S. A.

HUNTER

ET

AL

The conversion of energy flow of the storage element ,u~dnJdt to the reaction affinity-flow pair A, d[,/dt is performed by a conservative (neither energy storage nor dissipation) element

Electrical energy storage is also coupled with energy flow due to material transport through the relationship

(Zi,W,)~ = Y, zi32 = Y&J, ( ) noting that charge is conserved in the chemical reaction. The addition of storage elements and energy converters to the energy flow graph dictated by the irreversible production of entropy is also shown in Fig. 2. (B)

BOND

GRAPH

REPRESENTATION

THERMODYNAMIC

OF

THE

MODEL

The energy flow diagram in Fig. 2 is the bond graph for the threecompartment model of glucose metabolism. The addition of standard bond graph notation enhances the topological graph by separating the laws of conservation, continuity and the constitutive equations. Detailed explanations of bond graph theory are found in the literature (Auslander, 1972; Oster et al., 19733; Karnopp & Rosenberg, 1975). The bond graph representation of glucose metabolism illustrated in Fig. 3 was obtained by inserting the appropriate bond graph elements into the energy flow graph of Fig. 2. In constructing the bond graph, the assumption was made that the electrostatic potential of each compartment remains constant, and the capacitive storage element was replaced by an effort source. The two-port, non-linear resistance R, is used to designate the p reaction, and a capacitor Ci symbolizes the storage element for species i. Intracellular sources Si are associated with the chemical potential of species i assumed to remain constant and the metabolic flow of material in the citric acid cycle and the glycerol phosphate shuttle as noted previously in this report. The term Sf is the source associated with extracellular species i concentration. Bond graph signals not previously defined include the rate of intercompartment transfer to species i (JF), storage flow of species i (Ji), chemical reaction velocity (Jf) of the p reaction, the material flow of species i resulting from the p reaction (Jyp), the time rate of change of electrical charge (Jr), the entropy production of the p reaction (S,,), and the affinity for species i (At).

cw

/ C.,

U‘,,

JU

-

SW

1

0

TJ6,

At,. J. I If:,

A,,,

1J’ ,-

TF:-I

BOND (C)

GRAPH

DERIVATION

SIMULATION

OF

OF CONSTITUTIVE

STANDARD

CHEMICAL

METABOLISM EQUATIONS

509 AND

THE

POTENTIAL

The bond graph is based upon a theoretical analysis of the energetic interactions occurring within the system and does not involve considerations as to the nature of the constitutive equations relating the effort and flow variables of the particular bond graph elements.

(D)

CONSTITUTIVE

EQUATIONS

FOR

THE STORAGE

ELEMENTS

The capacitor corresponding to species i was defined as the partial derivative of the integrated flow variable with respect to the effort variable

Assuming that species i behaves ideally (metabolites are present at low concentrations), the chemical potential is related to the concentration ci by the relationship pi = &(I’,

T)+RT

In ci.

The term py(P, T) is the standard chemical potential of species i at a particular pressure P and temperature T The universal gas constant is represented by R. Using this relationship between the chemical potential and the concentration, the capacitance Ci was determined to be (Oster, Perelson & Katchalsky, 1973)

(E)

CONSTITUTIVE

EQUATIONS

FOR

THE RESISTANCE

ELEMENTS

Derivation of the non-linear resistance for each individual biochemical reaction involved defining the mechanism of the reaction and applying principles of equilibrium and/or steady-state enzyme kinetics to obtain a relationship between the reaction velocity and the substrate and product affinities. Numerical values for dissociation constants (K,), Michaelis coefficients (K,i) for species i, and forward and reverse rate constants (K,, K-& were reported in the literature or calculated from experimental results reported in the literature (Hunter, 1975). The resistance associated with the transport of lactic acid was based upon pure diffusion.

S. A.

510 (F)

CALCULATION

OF THE

HUNTER

ET

STANDARD

AL. CHEMICAL

POTENTIAL

The energy dissipated by a chemical reaction proceeding in the steadystate at constant temperature is the product of the reaction velocity and chemical affinity gradient. A numerical value for the affinity gradient requires an evaluation of the chemical potential, an intensive property of the reaction system, for both the reactants and products. Assuming ideal elements, the chemical potential (pi) is related to the concentration of species i via /.~i = ~(y(f’. T) + RT In Ci. The standard chemical potential (pi(P, T)) is independent of the concentration of species i and is identical to the standard free energy of formation ($) of the element i. Numerical values for the standard chemical potential were obtained by the simultaneous solution of linear, algebraic equations derived from the relationship where (AC”‘) represents the standard free energy of the chemical reaction. This analysis resulted in the formulation of 26 equations consisting of 35 variables, 26 of which were unknown. The results of the equilibrium thermodynamic analysis are summarized in Table 1. TABLE

Results of the equilibrium

1

thermodynamic study

(pi@‘, T) = Kcal mot p& phGATP phGADp p& p& /~&a &jHAp /i&P pkAD p;, pbpGA phADH &pGa jt+ ;;;; & /I&.

= = = = = = = = = = = = = =

- 219.22 - 876.0 - 678.6 - 420.6 -419.9 - 620.6 - 308.3 - 306.5 - 722.6 - 261.5 - 552.6 - 736.6 - 359.9 - 358.9

: 1 ;‘I’::: = - 133.9 = - 565.4

‘) piMP /A& $5,”

1 s,“,‘::;

/I& &s /& pkMN /I&, /&,, p&,, p&a /.l;Gpp, “2

= - 368.8 = - 419.4

= = = = = = = = =

-

695.8 761.3 108.99 419.8 56,69 158,3 300.2 328.5 574.1

I 1;‘;;: /I;; = - 2190 /I; = - 297.0

BOND

(G)

DIGITAL

GRAPH

COMPUTER

SIMULATION

SIMULATION

OF

METABOLISM

OF

GLUCOSE

511 METABOLISM

The response of the metabolic system to a change in plasma epinephrine concentration was simulated on a Digital PDP-10 system. The dynamic equations of the system resulting from the storage elements-were transformed to discrete equations and integrated with a simple step technique. The sampling period used was 10 ms. The simulated output included the Iconcentration and chemical potential of each metabolite, the velocity and entropy production of each chemical reaction, and the total entropy production of the system. 3. Results of the Computer Simulation The transient response of the model to a change in plasma epinephrine concentration from an initial value of 1.09 x lo-’ mol 1-l to a value of 2.73 x lo-’ mol 1-l was simulated. The response in terms of glycogen, glucose-l-phosphate, UDP glucose, glucose-6-phosphate and glucose concentrations are presented in Fig. 4 ; the transient response of the cyclic

0 GIP x IO3

0

40

80 Tme

120

(s)

FK;. 4. Transient response of the theoretical model to a change in plasma epinephrine concentration. The initial epinephrine concentration was 109 nM and the final concentration was 273 nM. CAMP, cyclic adenosine monophosphate; GLY, glycogen; G6P, glucose& phosphate; UDPG, uridine diphosphoglucose;GLU, glucose;GIP, glucose-l-phosphate.

3’,5’ adenosine monophosphate-dependent protein kinase control system is shown in Fig. 5. Other model metabolites displayed only a small change from the initial values during the 2 min simulation. The phosphoglucomutase reaction changed direction from glycogen synthesis to glycogen degradation between 5 and 10 s after the change in epinephrine concentration occurred. The phosphorylase a reaction velocity, initially at

512

S. A.

HUNTER

ET

AL.

FIG. 5. Simulated response of the CAMP-dependent control system to a change in epinephrine concentration. Epinephrine concentrations as in Fig. 4. Phos b, phosphorylase b: Phos a, phosphorylase a; APhos K, active phosphorylase kinase; IPhos K, inactive phosphorylase kinase; Phos ratio, (Phos a)/(Phos bl; Phos K ratio, (APhos K)/(IPhos K); GSI, glycogen synthase-independent.

3.2 x 10y8 mol mm’ g-’ of muscle, peaked to a value of 3.7 x lo- 6 mol min-’ gg’ of muscle 35 s after the disturbance occurred and then decreased to a value of 4.1 x 10-l mol min-’ g- ’ of muscle at the end of the 2 min simulation. With the exception of the hexokinase reaction velocity which decreased slightly, the velocities of reactions in the Embden-Meyerhof pathway increased in magnitude, phosphoglucoisomerase by a factor of 1.3. The entropy production resulting from the disturbance in plasma epinephrine concentration is illustrated in Fig. 6.

0

80

40 Time

(s)

FIG. 6. Calculated transient entropy production resulting from the disturbance in plasma epinephrine concentration given in Fig. 4.

BOND

GRAPH

SIMULATION

OF

METABOLISM

513

4. Evaluation of the Computer Simulation (A)

QUALITATIVE

EVALUATION

The simulated response to a step increase in plasma epinephrine concentration corresponds qualitatively to biochemical descriptions found in the literature (McGilvery, 1970). The epinephrine stimulation increased the adenyl cyclase reaction velocity by a factor of four resulting in a rapid buildup of cyclic 3’,5’ adenosine monophosphate. Ninety-seven percent of the cyclic 3’5’ adenosine monophosphate was bound to the regulatory unit of protein kinase, freeing the catalytic unit of protein kinase which catalyzed the conversion of glycogen synthase from the independent to the dependent form and phosphorylase kinase from the pH -6.8 inactive to active form. Active phosphorylase kinase, in turn, converted phosphorylase b to phosphorylase a, the final step in the cascade control system, causing the rapid breakdown of glycogen to glucose-l-phosphate. Increasing the concentration of glucose-l-phosphate changed the direction of the phosphoglucomutase reaction from glycogen synthesis to glycogen degradation. The rising concentration of glucose-6-phosphate inhibited the hexokinase reaction slightly, decreasing the reaction velocity by 2% at the end of the 2 min simulation. The inhibition of hexokinase was not overcome by the 17% increase in intracellular glucose concentration, a characteristic of mixed inhibitors. The increase in glucose-6-phosphate also resulted in an increase in the readily reversible phosphoglucoisomerase reaction and a buildup of fructose-6-phosphate. Other reaction velocities and metabolite concentrations in the EmbdenMeyerhof pathway remained fairly constant. The pathway was essentially shifting from a metabolism of glucose to a metabolism of glucose-6phosphate derived from glycogen degradation to maintain the necessary production of adenosine triphosphate. This reflected the basic idea that the pathway functions to maintain the level of adenosine triphosphate. The conversion of glycogen synthase from the independent to the dependent form proceeded at an extremely slow pace while the reaction velocity increased by a factor of two over that of the initial state. The reaction model is based upon converting the free, uncomplexed enzyme. The increase in uridine diphosphoglucose resulted in a decrease of the free glycogen synthase and a decreased amount of substrate for the cyclic 3’,5’ adenosine monophosphate-dependent protein kinase reaction. The fact that the reaction velocity increased by a factor of two in the face of a decreased substrate concentration indicates that the control system was activated. The entropy production response illustrated in Fig. 6 peaked 20 s after the input disturbance. The major component of this peak was the entropy

514

S. A.

HUNTER

ET

.4L

production due to the phosphorylase a reaction and represented the potential energy of glycogen dissipated during the degradation of glucose-lphosphate. The time rate of change of entropy production approached zero, and the entropy production became constant, the thermodynamic criterion for a steady-state. The entropy production response indicates a new steadystate resulted from the input disturbance. The response in Figs 4 and 5 also indicate the transient is dying out. The energetics of the protein kinase control system were not included in the entropy production analysis. (B)

QUANTITATIVE

ASSESSMENT

Stull & Mayer (1971) have performed in viva experiments on rabbit gracilis muscle and measured the response to an injection of 4 x 10m9 mol of the cathecholamine isoproterenol. Kinetic data of Severson er al. (1972) show that the activation of adenyl cyclase by isoproterenol is a good approximation of the activation by epinephrine. The magnitude of the epinephrine disturbance for the simulation was calculated from the average concentration of hormone in muscle measured by Stuli & Mayer (1971). Comparisons between some experimental results of Stull & Mayer (1971) and theoretical values obtained through the previously described simulation are shown in Figs 7 and 8. The transient increase in the cyclic 3’,5’ adenosine monophosphate increased by a factor of two after 60 s and remained constant as compared to a factor of approximately 3 f 05 in the experimental response. The simulated response had a phosphorylase activity

8

c

Phosphorylose T

3

E 8 2

Kmase

ratlo x IO2

Experimental data from WI and Mayer with the phosphorylase kmase rate x 06 x IO+’

4

0

80

40 Time

120

(s)

FIG. 7. Comparison of experimental (closed circles with standard deviation indices) and calculated values (open symbols) of phosphorylase kinase ratio and phosphorylase ratio following a step increase in isoproterenol or epinephrine, respectively. Experimental data from Stull 81 Mayer (1971); theoretical values as shown in Fig. 5.

BOND

GRAPH

SIMULATION

OF

METABOLISM

515

00000000000000000 Smhted response

0

40

80

120

Trne (sl

FIG. 8. Comparison of experimental (closed circles with standard deviation indices) and calculated values (open circles) of CAMP following hormonal stimulation. Experimental data from Stull & Mayer (1971); hormones as in Fig. 7.

ratio of O-25 after 90 s in comparison to a value of approximately 0.2 + 0.1 in the experimental response after correction for the non-cyclic 3’,5’ adenosine monophosphate value. The change in the phosphorylase kinase activity ratio was equal to a value of 0.02 whereas the change in the activity ratio of the experimental response was OG9 &- 0.01. Other pieces of information have been published which tend to validate the simulated response. Drummond et al. (1969) injected 5 c(g of epinephrine per kilogram of rat intracardially and measured the response of rat gastrocnemius muscle after an elapsed time of 40 s. Drummond et al. (1969) reported the following changes : (1) increase in the concentration of cyclic 3’,5’ adenosine monophosphate from 0.3 to 1.6 pmol per kg of muscle (5.08 x lo-’ M to 2.71 x 1O-6 M); (2) increase in phosphorylase - AMP/+AMP ratio from 0.25 to 0.62 ; (3) increase in phosphorylase kinase 6.8/8.2 ratio of 0.03 to 0.13. Posner et al. (1965) injected 10 pg of epinephrine per kg of rat intracardially and measured the response of rat gastrocnemius muscle after an elapsed time of 60 s. Posner et al. (1965) reported the following results: (1) increase in cyclic 3’,5’ adenosine monophosphate from 069 + 0.5 to 1.82 + 0.75 pmol per kg muscle using one approach and from 0.66 & 0.05 to 1.2 + 0.42 pmol per kg muscle using a second method ; (2) increase in phosphorylase - AMP/+ AMP ratio from O-1 f 0.04 to 048 + 0.13; (3) increase in phosphorylase kinase 6.8/8.2 ratio from 004 f 0.04 to 0.25 + 0.12. 19

516

S. A.

HUNTER

ET

AL

In experiments with a perfused rat heart, Cheung & Williamson (1965) obtained data which indicates that cyclic 3’5’ adenosine monophosphate is not an important activator of phosphofructokinase in vim. This lack of regulation was also present in the simulated response. Helmreich & Cori (1965) reported the results of experiments performed on isolated, intact frog sartorius muscle subjected to a 5 x lo-’ M epinephrine disturbance. Their experiments showed that a large increase in glucose-6-phosphate and fructose-6-phosphate was not necessarily accompanied by a corresponding increase in phosphofructokinase activity. The basal rate of lactate formation in those experiments was 0.01 umol of lactate per ml per min. During the first 30 min of epinephrine action, when the level of hexose monophosphates had risen nearly lo-fold and the phosphorylase a level had presumably increased to 40% of the total, there was no increase in lactate formation. The observation thus tends to support the computer simulated response of a constant metabolic throughput from phosphofructokinase through the rest of the Embden-Meyerhof pathway. 5. Summary A dynamic simulation of skeletal muscle glucose metabolism has been presented. The mathematical formulation consisted of a thermodynamic analysis of irreversible processes in conjunction with the dynamics of energy storage by reversible elements. The bond graph topology, isomorphic with respect to the system physiology, was based on the biochemical reactions of the Embden-Meyerhof pathway and the protein kinase controlled synthesis/degradation of glycogen stores. The constitutive equation for each biochemical reaction was developed from experiments reported in the literature. These non-linear algebraic equations were derived using equilibrium and steady-state enzyme kinetics. Constitutive equations for the capacitive storage elements were based upon the assumption of ideal metabolites. The coefficients used in the model construction were developed on physical principles as opposed to coefficients computed from input-output data in a “black box” analysis. The modeling process led to the evaluation of unknown constants such as the standard chemical potentials of species in the biochemical pathways under investigation and the chemical reaction dissociation constants for several enzymatic reactions. The modeling process, in addition, resulted in questions such as those concerning the location and effective concentration of enzymes. Essentially, the modeling process forced an evaluation of the current-state-of-the-art and an extension of existing theory and/or information.

BOND

GRAPH

SIMULATION

OF

517

METABOLISM

The application of thermodynamics and bond graph theory to metabolism resulted in a model which simulated the response of muscle to a plasma change in epinephrine concentration. The response was in qualitative and quantitative agreement with metabolic changes.

glucose skeletal model known

REFERENCES ALLEN, R. AUSLANDEQ.,

R. ( 1978). Comp. Prog. in Biomed. 8, 149. D. M., OUTER, C. F., PERELSON,

A. &

CLIFFORD.

C. (1972).

Trans. Am. Sot. Mech.

Eng., J. Dynamic Sysr., Meas. Cont. 94, 239. CHEUNG, W. Y. & WILLIAII(SON, J. R. (1965). Nature 208, 970. DRUMMOND, G. I., HARW~~D, J. P. & POWELL, C. A. (1969). J. biol. Chem. 244,4235. GEBBEN, V. (1977). Trans. Am. Sot. Mech. Eng., J. Dynamic Syst., Meas. Cont. 99, 143. HELMREICH, C. & CORI, C. F. (1965). Adv. Enz. Reg. 3,91. HUNTER, S A. (1975). Ph.D. Thesis, Univ. Microfilms, Ann Arbor, Michigan. INCROP~A, F. P. (1975). J. Environ. Qual. 4, 440. KARNOPP, D. C. & ROSENBERG, R. C. (1975). System Dynamics: A Unified Approach. New York: Wiley-Interscience. LAKE, J. V (1967). Aust. J. Biol. Sci. 11, 487, 494. LOMMEN, P. W., SCHWINTZER, C. R., YOCUM, C. S. & GATES, D. M. (1971). Plnnra 98, 195. MCGILVERY. R. W. (1970). Biochemisrry: A Functional Approach. Philadelphia: W. B. Saunders co. OSTIZR, G. F., PERELSON, A. S. & KATCHALSKY, A. (1973). Q. Rev. Biophys. 6, I. POSNER, J. B., STERN, R. & KREBS, E. G. (1965). J. biol. Chem. 240, 982. PRICoCiINE, I. (1967). lntroducrion to Thermodynamics I$ Irreversible Processes. New York : Interscience Publ. RAPOHIRT, D. L., HEINRICH, R. & RAPOPORT, S. M. (1976). Biochem. J. 154,449. SEVERSON. D. L.. DRUMMOND. G. I. & SULAKHE. P. V. (1972). J. biol. Chem. 247. 2949. QIJLL, J.‘T. & GAYER, S. E. i1971). J. biol. Chem. 246; 5716. VON BERTALANFFY, L. (1962). Trends in General Systems Theory, p. 24. New York: John Wiley & Sons.