Energetics of Glucose Metabolism: A Phenomenological Approach to Metabolic Network Modeling

A new formalism to describe metabolic fluxes as well as membrane transport processes was developed. The new flux equations are comparable to other phenomenological laws. Michaelis-Menten like expressions, as well as flux equations of nonequilibrium thermodynamics, can be regarded as special cases of these new equations. For metabolic network modeling, variable conductances and driving forces are required to enable pathway control and to allow a rapid response to perturbations. When applied to oxidative phosphorylation, results of simulations show that whole oxidative phosphorylation cannot be described as a two-flux-system according to nonequilibrium thermodynamics, although all coupled reactions per se fulfill the equations of this theory. Simulations show that activation of ATP-coupled load reactions plus glucose oxidation is brought about by an increase of only two different conductances: a [Ca2+] dependent increase of cytosolic load conductances, and an increase of phosphofructokinase conductance by [AMP], which in turn becomes increased through [ADP] generation by those load reactions. In ventricular myocytes, this feedback mechanism is sufficient to increase cellular power output and O2 consumption several fold, without any appreciable impairment of energetic parameters. Glucose oxidation proceeds near maximal power output, since transformed input and output conductances are nearly equal, yielding an efficiency of about 0.5. This conductance matching is fulfilled also by glucose oxidation of β-cells. But, as a price for the metabolic mechanism of glucose recognition, β-cells have only a limited capability to increase their power output.


Introduction
In recent years mathematical modeling of metabolic pathways became increasingly important as a helpful tool for a deeper understanding of biochemical reaction networks in cellular metabolism and energetics. Besides older kinetic approaches, several new models were developed, whose flux equations contain driving forces, primarily to eliminate possible inconsistencies with the Second Law of Thermodynamics [1][2][3][4][5][6][7][8][9]. In two recently published articles [6,10], a new flux equation was introduced, which allows both, the description of kinetic flux behavior similar to Michaelis-Menten type equations [11], as well as the treatment of energetic relations known from methods of

Flux Equations
In previous articles (Diederichs,[6,10]) a new phenomenological type of flux equation to describe not only membrane transport reactions but also enzyme-catalyzed reactions as they occur in the living cell, was introduced. The new formalism combines some characteristics known from enzyme-catalyzed reactions with energetic relations of NET. They are comparable to several other linear relationships between corresponding flows and forces such as electric current flow, heat flow, or diffusional flow. The main difference to these latter phenomenological laws and to NET is brought about by introducing a variable instead of a constant conductance into flux equations. Conductances may be varied by substrate and ion concentrations as is known from enzyme kinetics [11], leading to the following expression:  is the dissipation function per unit volume, it is given in units of power dissipated per unit volume. It is identical to the heat dissipated by the irreversible reaction per unit time and volume. Analogous to electrical power P = I  U and current I = L  U, the expression of flux J then is given by Equation (1a). This latter formulation, which is used throughout this paper, may be advantageous in view of its similarity to relations known from techniques such as electric circuits and heat flow. It should be noticed, however, that A is not identical to a conjugated force X, which according to Equation (2a) is defined by X = A/T. In the following, the relationship to a Michaelis-Menten type reaction is shown. The flux through an enzyme-catalyzed reaction may be expressed as       The S (S = substrate or effector) containing term of Equation (3a) and (3b) represents a dimensionless activation/inhibiting factor, which can be regarded as a concentration-dependent probability, by which the maximal conductance max Enz L of an enzyme-catalyzed reaction may be reached.
If several activating and/or inhibiting factors were involved with the catalytic process, then, according to probability calculus, these factors must appear as products in the respective flux Equation [6].
On the other hand, Equation (3a) can be regarded as a special case of flux equation derived from the more fundamental rate Equations (1a) or (3b). In Equation (3a) the variable term max . The former contain variable conductances, but constant affinities, whereas the opposite is true for equations of NET; these contain constant conductances, but variable affinities.
The great advantage of flux equations containing both variable conductances as well as variable affinities, is that they are remarkably well-suited for metabolic network modeling. So flux control by activating or inhibiting effectors is made possible as it is known from enzyme kinetics [11]. The inclusion of driving forces has the advantage that energetic parameters like power output or efficiency of coupled processes as they occur in metabolic pathways can be calculated.
As was shown in reference [6], the new formalism is applicable to reactions under "far from" as well as under "near" equilibrium conditions. Far from equilibrium, the flux equation is dominated by activation factors. This is brought about by the logarithmic nature of A, which remains fairly constant in large logarithmic arguments. Under these conditions, the flux equation approaches Michaelis-Menten kinetics [6] and, therefore may be suitable to describe far from equilibrium reactions, which is demonstrated by present simulations (most reactions are far from equilibrium). Theoretically, it is possible to reverse such an irreversible flux, but only at extreme substrate and/or product concentrations which do not occur in living cells.
Under near equilibrium conditions, the flux equation becomes dominated by A, because the logarithmic function then approaches zero. In this region the slope of this function becomes markedly increased. Under these conditions, the flux equation approaches the NET formalism. This approach is all the more complete, the more the K M value of this reaction approaches zero.
In Figure 1

Constant versus Variable Conductances
To demonstrate the different behavior between pathway fluxes composed of constant and variable L's respectively, a metabolic pathway composed of n in series enzyme-catalyzed reactions may be regarded. If, e.g., L i is increased by a factor f at constant L's, the pathway flux J expectedly would be increased by a factor f′ < f. Then, according to Ohm's law, the fractional increase of J is equal to the fractional increases of all n − 1 affinities at the n − 1 remaining conductances. Only the fractional change of A i at L i (fractional increase of L i is imposed and given by f − 1) must be different from all others. It is given by A/A 0 = (J/J 0 − (f − 1))/f. The above conditions of constant conductances are equivalent to a simple electric network, but are probably not compatible with a metabolic pathway. Since, under these latter conditions, variable enzyme activities determine conductances and thus might be responsible for a quite different behavior. The n − 1 fractional changes of affinities are not equal as is to be expected for constant conductances. This fact is associated with a remarkably different quality. Pathway fluxes containing variable conductances, after a given perturbation, can reach a new steady state much more rapidly than those having constant conductances. This is verified with a simple simulation (SIM GLY , see App. A14, the maximal conductance of the first reaction step, max GK L , is increased by a factor f = 2.0    In this context, it seems worth mentioning that Michaelis-Menten type formulations of pathway flux (const. A's to yield max max .

Enz
const Enz , under the same conditions are unable to produce a new steady state. Obviously, in addition to variable conductances, the inclusion of variable driving forces in pathway fluxes is indispensible for the toleration of non-infinitesimal changes of conductances as they occur in living cells, for instance during activation of glucose metabolism, to increase power output.
With the new formalism, however, it is possible to describe metabolic networks as was shown recently with the known phenomena of stimulus-secretion coupling of β-cells [6,10]. In simulations the increase of only one single parameter, the glucose concentration, is necessary and sufficient to induce β-cell activation, as is known from experiments.

Phenomenological Relations
In a metabolic pathway, which may be composed of a source and a sink and several in series enzyme-catalyzed reactions at steady state, the pathway flux can be expressed as This follows from  = J i A i , and at steady state  = JA i = JA ov . Furthermore, due to the similarity of the flux equation to Ohm's law, in series conductances generate an equivalent conductance L ov . A constant overall affinity of the pathway A ov = A i is then determined by the constant concentrations of metabolites of the source and sink, respectively. So, it is possible to contract several in series reactions of a given reaction sequence to one single step. In cellular metabolism this principle is used for metabolic channeling [32,33], which is brought about by aggregating several different enzyme molecules of a sequence to one catalytic complex. Because respective intermediates do not appear in solution, this may cause an appreciable increase of the overall conductance of the catalytic complex.
Since chemical and electro-chemical potentials are potential functions, their line integral over a closed path must vanish, i.e., Kirchhoff's loop rule is applicable also to metabolic networks [2]. Kirchhoff's second rule, the junction rule, however is not always valid. For instance, the aldolase reaction splits the first part of GLY into two equal fluxes, which merge again to the second part of GLY. The flux of this latter pathway is not equal to that of the first part, as is to be expected for an electric network (conservation of charge), but is increased by a factor of two. If, however a given flux is split up into several parallel fluxes, the junction rule is obeyed [2]. In any case, the law of mass conservation has to be fulfilled.
This latter fact has to be taken into account, when A ov and L ov are to be determined. Equation (4) is applicable only if a constant flux of given magnitude through the entire pathway is realized. With regard to GLY this would mean that for instance the flux of the second part (J Ga , A3) has to be transformed to yield the same magnitude as the flux of the first part. Since a given  can be expressed by any other reliable product of J and A, and consequently also by a product containing a J of desired magnitude, e.g., J GK (which is a measure of glucose utilization), an altered affinity has to be determined in such a way that  remains unchanged. This transformed affinity is given by and the associated transformed conductance by Such transformed values have to be used, when the steady state pathway flux of a metabolic network containing one or several fluxes of different magnitude is to be determined. To develop a simulation of a metabolic network like coupled glucose oxidation, it is not necessary to know all parameters exactly. If glucose utilization and/or oxygen consumption are known from experiments, conductances can be adjusted for a maximal power output. Affinities of single reactions often are given or may be estimated from Gibbs energies of formation [Alb]. With known stoichiometric coefficients the simulation must yield correct results also with respect to thermodynamic constraints, even if the steepness of flux equations is not precisely met. This means that the use of adjusted constants like K M values and other modeling parameters may suffice to fulfill reliability requirements.

Flux Control Coefficients
In Metabolic Control Analysis (MCA) several dimensionless coefficients are defined, to quantify metabolic control at steady state [34]. One of these coefficients is the flux control coefficient J Ex C , which is defined as (E x = one particular enzyme concentration of a metabolic pathway of several reaction steps in series; J x = change of pathway flux produced by a perturbation E x /E x = f − 1). According to MCA, the sum of flux control coefficients ( first part of GLY (J GK , J Ald , and J Ti ) are equal and a new steady state pathway flux is produced. All new fluxes are in fact increased, but the relative increase is much less than f − 1 (0.0694 compared to 1.0).
For conditions of constant conductances, the same simulation as described above but with constant L's is used. For these conditions, the sum of J Enz C can be obtained analytically. It is derived from  For variable L's, the above relation cannot be solved analytically. The output values of the simulation, however, show that respective points follow roughly the analytical function. At small perturbations ( 1.0 f  ) the summation theorem seems to be sufficiently approached, as is shown above for f = 1.01. For larger changes, however, as they normally occur in cellular metabolism, deviations from 1.0 have to be expected ( Figure 3). In a more complex pathway like whole coupled glucose oxidation, coupled reactions and, in addition, load reactions ('load' usually designates the output affinity A 1 of a cytosolic reaction coupled to ATP c splitting) are involved, both of which may impair the constancy of A ov . As a result, fractional changes of the pathway flux are also influenced by a changed A ov . This is given by According to Equations (5a) and (5b), transformed values have to be used, which are designated tot A  and tot L  , respectively, to describe more complex networks. Such influences become most pronounced, when uncoupled reactions are involved (see Section 2.3.1).

Energetic Coupling in Metabolism
In metabolism, the chemical free energy of one mole of hydrogen containing substrates like glucose is converted ultimately into transformed Gibbs energy of reaction of cytosolic ATP hydrolysis ( r G′(ATPc) = −A ATPc ) plus molar heat changes. This energy conversion is brought about by reactioncoupling and is widely described by NET [12][13][14][15]. The following known equations are usually used to describe coupling of a two-flux-system: J 1 and J 2 designate the output and input fluxes (related here to d/dtV), respectively, X 1 and X 2 represent the output and input forces; here they are identical to driving forces or affinities A 1 and A 2 , respectively. L 11 and L 22 are termed straight coefficients, and L 12 is termed coupling coefficient. The degree of coupling is defined as 12 11 22 q L L L   (8c) the phenomenological stoichiometry as 11 22 To show more transparently, how coupled and uncoupled fluxes work together, the above equations are given in an equivalent formulation for enzyme-catalyzed reactions and associated free energy conversions occurring in cellular metabolism [10]: L c represents the coupling conductance (= L 12 ), whereas  1 and  2 are related to the leak conductances L L1 and L L2 by  1 = L L1 / L c , and by  2 = L L2 / L c .
The output affinity A 1 is usually opposed to the input affinity A 2 , i.e., A 1 is negative, their sum however must be positive. At non-zero leak conductances, leak fluxes are produced in the direction of their respective driving forces. They are given by J L1 = L L1 A 1 (negative), and J L2 = L L2 A 2 . Comparison of Equations (8a) and (8b) with Equations (8f) and (8g), yields L 11 = L c + L L1 , and L 22 = L c + L L2 .
At q = 1.0 ( 1 and  2 are both zero, see Equation (8j)), J 1 and J 2 are identical and are given by . The resultant driving force of in series affinities is given by When  1 and/or  2 > 0 (q < 1.0), it can be taken from Equations (8f) and (8g) that now fluxes must be different, and that J 1 becomes decreased by the negative leak flux J L1 , whereas J 2 is increased by J L2 . Obviously, uncoupling is always brought about when such additional leak fluxes arise.
The degree of coupling in terms of 's can be derived from Equations (8f) and (8g): The root has to be drawn, because from both flux ratios, J coup1 /( J coup1 + J L1 ) and J coup2 /( J coup2 + J L2 ), associated with affinities A 1 and A 2 , respectively, only one resultant ratio can be associated with the flow through the coupling conductance L c . Thus, q represents the geometric mean of both ratios involved.
The efficiency of free energy conversion is given by with x = (A 1 /A 2 )Z (reduced force ratio). Maximal efficiency is given by (8n) and efficiency at maximal power output by 2 max 2

Coupling in Oxidative Phosphorylation
During oxidative glucose metabolism hydrogen of glucose is transferred to oxidized hydrogen carriers like NAD ox , which in turn transfer their hydrogen to the respiratory chain of the mitochondrial inner membrane, where it reacts with oxygen to form H 2 O and reoxidized carrier. This reoxidation of redox carriers (NAD red and FAD red , respectively) at the inner side of the mitochondrial inner membrane is coupled to the transport of protons from the mitochondrial matrix space to the outer side of the inner membrane. At the interface between the outer aspect of the inner membrane and the bulk solution of the intermembrane space (solution between inner and outer mitochondrial membranes) these protons may exchange with other cations like K + ions of the intermembrane space, so that an electrochemical potential difference of protons over the inner membrane, H   , is produced by coupled proton transport. It is composed of a chemical potential difference of protons over the inner membrane,  H , plus an electrical potential difference, the membrane potential  m  of the inner membrane, For simulations presented in this and foregoing articles [6,10] it was assumed that the pH of the intermembrane space is identical to the cytosolic pH c and that outwards pumped and inward flowing protons only affect matrix pH m , but have no effect on pH c (e.g., J NA or J SY , A7 or A9). The same holds for proton symport and exchange reactions at the inner membrane: only pH m can be changed, pH c remains unaffected (e.g., J Pi and J HC , see [2]). In contrast,  m  is changed by both, in -and outgoing electrogenic membrane fluxes like J CU (see [2]) or J SY .
Such behavior might be realized for real mitochondria in situ by a further compartmentation of the intermembrane space. From recent morphological studies [35][36][37][38], it is known that cristae membranes form flat sacs with narrow openings to the intermembrane space. If proton cycling (outwards pumping and inwards flow) would preferentially occur at these structures, then primarily the luminal pH of cristae would be affected. This means that H   of the inner membrane might be transformed into a potential difference with a more pronounced (more negative) chemical component ( H ) and an appropriately less pronounced (less negative) electrical component ( m ) at cristae membranes than at the rest of the inner membrane. A low buffering power of these regions would further increase this effect. Strauss et al. [36] could demonstrate that ATP synthase complexes are located preferentially at the rims of cristae sacs where they may be involved with shaping the rim structure. These morphological facts are consistent with the idea that pH c under physiological conditions may be affected only to a negligible degree by proton fluxes across the inner membrane. So, an amount of H n    (J) of energy is dissipated, when n moles of protons move without coupling across the inner membrane from the bulk phase of the intermembrane space to the bulk phase of the matrix. Redox reactions (J NA and J FA , A7 and A8) are coupled to proton transport with a stoichiometry of 10.0 and 6.0 mol H + per mol of NAD red or FAD red oxidized [18], respectively. Thus, 6.0 mol O 2 and 12.0 mol reducing equivalents are consumed per mol of oxidized glucose. These chemical constraints are fulfilled by all simulations. But not all protons transported by J NA and J FA enter oxidative phosphorylation (OP). Some leak back into the matrix (J PL ), some are needed for Ca/H exchange (J HCE ), and some for the malate-aspartate shuttle (J MS ). In addition, in real mitochondria H + influx may be required for transhydrogenation. This membrane reaction, however, is not addressed in present simulations. The remaining proton flux is used by ATP synthase (J SY ) and ATP/ADP exchange (J AE ) plus H/Pi symport (J Pi ). At steady state, the stoichiometry of coupled J SY is assumed to be 3.0 mol H + per mol of matrix ATP (ATP m ) produced [18]. ATP m is transferred into the cytosol and ADP c plus Pi are delivered to the matrix by two further in series reaction steps (via J AE and J Pi , respectively). Both reactions are coupled to inward transport of 1. were included, this value would be slightly higher).
When ATP splitting by load reactions is added, the whole process may be considered as two interconnected cycles: a proton cycle coup H J coupled to an ATP cycle coup ATP J through the ATP synthase reaction, ATP/ADP exchange, and H/Pi symport. At steady state (all individual fluxes of a given cycle are of same magnitude), both cyclic flows produce no dissipation (or entropy), because their respective sums of affinities taken over a closed path in both cases must be equal to zero.
In the above derivation proton fluxes from proton sources to proton sinks at the outer and inner sides of the mitochondrial inner membrane are not included in calculations. Under real conditions, however, in order to allow proton cycling, such fluxes must be present, and thus, also a certain amount of dissipation of free energy. In addition, during ATP cycling free energy is dissipated by diffusional fluxes of ATP, ADP, and P i species (see Section 2.4). It can be concluded, therefore, that under real conditions OP must be always slightly uncoupled, even if all respective inner membrane reactions of OP would proceed at q = 1.0. The

Indirect Uncoupling
Equation (10) cannot give any information about the mechanisms of coupling and uncoupling, respectively. However, from their derivation [39] it follows that uncoupling leak fluxes must flow through the coupling device, i.e., through the respective multi-enzyme complex involved, because they are added to input and output fluxes, respectively. In contrast to this direct uncoupling an indirect uncoupling may occur through additional parallel leak fluxes like the proton leak flux J PL at the inner mitochondrial membrane, or the flux of Na + and K + ions through sodium and potassium channels or other Na + and/or K + -permeable channels of the cell membrane. The relevant coupled reaction for these latter uncoupling fluxes is given by the reaction sequence of the Na-K pump. Also Ca 2+ pumps must be mentioned in the context of such pump and leak fluxes. Remarkably, dissipation of free energy is used further by these systems for signal transduction.
In principle, to avoid static head production, coupled reactions are controlled in addition to parallel leak fluxes by their own substrates. For instance, Na-K pumps are kinetically controlled by [Na + ] c , and Ca 2+ pumps by [Ca 2+ ] c , so that pump fluxes slow down when substrate concentrations become sufficiently lowered.

Stucki's Conductance Matching
Stucki suggested that energy coupling through OP may be described as a two-flux-system in the same way as for a single coupled reaction [40]. He identified the oxygen utilizing redox reaction as the input and the ATP c yielding reaction as the output flux. Including an additional attached ATP c consuming load reaction, he could show that the dissipation function (entropy production in his article) of such a system possesses a minimum at a reduced force ratio given by 33 and L 11 are Stucki's load and input conductances, respectively). Equating this value with the x coordinate of maximal efficiency, , shows that a minimal dissipation and a maximal efficiency would coincide, if the ratio of the load (L 33 ) and the input (L 11 ) conductance of whole OP fulfilled the condition 2 33 11 1 L L q   . This relation was termed conductance matching by Stucki [40]. It is derived for uncoupled systems only, because under totally coupled conditions the above relationship requires that L 33 vanishes. This, however, would generate a static head state, which is not compatible with conditions of metabolizing cells.
It seems already worth mentioning that all simulations presented here and those published recently [6,10] do not adjust to such a value. Moreover, they all allow total coupling without producing a static head.
To concluded that for a given q ov (overall value of q for whole OP), X 1 may be the only variable of this two-flux-system, and that all other parameters, i.e., X 2 , L 11 , and L 12 are constant and can be derived from experimental results [40].
The simulation demonstrates, however, that linearity under such conditions can be produced with variable conductances and variable forces as well. It can be taken from Figures 4 and 5 (panels A and B) that the fluxes J 11 and J 12 , whose sum yields J 1 , are slightly curved with opposite curvature, but that their sum yields a straight line. The same holds for J 2 . The slight deviation from linearity shows that some results of the simulation only approximately obey flux equations of NET. Under coupled conditions, fluxes J 1 and J 2 are equal and linear (not shown).   Figures 4d and 5d demonstrate that uncoupling of individual reactions of OP fulfill the analytical efficiency equation (10). Panels D-F of Figure 4 show that for this individual two-flux-system (A NA is far from equilibrium) efficiency is close to the maximum, but that dissipation is not minimal, whether at the reduced force ratio x max (where efficiency is maximal), nor for that obtained from the simulation. Power output of this reaction is 67.5% at q n = 0.952. In contrast, panels D-F of Figure 5 show that efficiency of this near equilibrium reaction (A SY is near zero) is not maximal. Dissipation is close to a minimum for the x value obtained from the simulation, but not for the  max related x. As can be expected, power output for this near equilibrium reaction is much lower, only 11.7% at q s = 0.971.
The It must be zero, because the output force is zero under all conditions. A conductance matching, as demanded by Stucki, therefore does not exist in simulations presented here, and, since the reaction sequence of ATP formation and ATP consumption must be present as an indispensible part of energy metabolism, this mechanism may also not be realized in a living cell.   Furthermore, from a teleological perspective, Stucki's conductance matching seems highly unlikely to occur in energy metabolism of living cells, since such behavior would lead to the paradoxical situation that most favorable conditions of energy conversion -namely power output without losses of free energy by uncoupling -must result in cell death by the inevitable emergence of a static head.
When A , respectively) must vanish to describe the involved reactions as a two-flux-system. Although this is possible, unreliable results are produced, which often show negative ov  's and a q ov > 1.0, respectively. This might be caused by the cyclic mode of operation of the coupling system. It is concluded, therefore, that it is impossible to describe OP as a two-flux-system according to NET, although individual reactions of this complex coupling system may fulfill those equations.

Activation of Power Output in Glucose Metabolism
In the following, how the pathway of glucose oxidation may react on increasing demands of power output including uncoupling will be demonstrated. Two very different cell types-the ventricular myocyte with a high potential of power output, and the pancreatic β-cell with a relatively low power output-will be compared.

Ventricular Myocytes
Hexokinase of heart muscle has a very low K M value compared to that of glucokinase of β-cells. It is this kinetic difference of the first glycolytic reaction step, which is necessary to allow high pathway fluxes of glucose oxidation in muscle cells. In addition, the low K M prevents that the intracellular glucose concentration (  (7)).
Therefore, if a J Enz C  >> 1.0 were found with a pathway containing coupled reactions, this may indicate a significant uncoupling. The above results of control analysis show that increasing only two conductances may be sufficient to appreciably increase pathway flux and thus, also power output. A simultaneous increase of all maximal conductances is not necessary; the variability of conductances obviously is able to produce the demanded increase of whole pathway flux. Moreover, also fluxes through J AE and J Pi must not be specifically activated. Ca 2+ activation of pyruvate dehydrogenase (PDH), and dehydrogenases of the citric acid cycle (fluxes J PDH and J CAC , respectively) is virtually without effect on power output and O 2 consumption in present simulations, because respective conductances ( max In real VMs, these conductances might be appreciably lower, so that a [Ca 2+ ] m activation becomes indispensable at an increasing power output [23,24,[41][42][43]. Also ATP synthase of VMs is known to be activated by [Ca 2+ ] m [44]. In present simulations, activation of ATP synthase (J SY ) by [Ca 2+ ] m cannot increase power output, but reduces dissipation of free energy of this reaction. Beside these latter effects of [Ca 2+ ] m , mitochondrial Ca 2+ fluxes may be needed to synchronize free energy flows of several mitochondria during contraction cycles. Otherwise an optimal interplay of all involved reactions might not be guaranteed.
In the present simulation as well as in real VMs, important load reactions like muscular contraction or Ca 2+ transport by SERCA pumps are activated by an increase of [Ca 2+ ] c . This in turn leads to an increase of ATP c consumption with a concomitant elevation of [ADP] c and [Pi] c . Because the adenylate kinase (AK) reaction is near equilibrium, this means that also [AMP] c must increase and PFK becomes activated, with the result that all pathway fluxes associated with glucose oxidation become activated too. So activation of ATP c utilization is followed by an appropriately adjusted ATP c delivery. In this way load reactions not only deliver useful work, but concomitantly generate [AMP] c , which quasi as a third messenger activates ATP c production, so that a drastic decrease of [ATP] c and associated A ATPc may be prevented, even under conditions of a markedly increased power output. Present simulations demonstrate that the increase of only [Ca 2+ ] c dependent load conductances and [AMP] c dependent L PFK is needed, to switch from low to a high power output. O 2 consumption of working hearts shows linear dependence on mechanical power output over a wide range of delivered power [23][24][25]. Such behavior can be simulated, supposing that pressure development of whole hearts is mechanically produced by the contractile machinery of single VMs, i.e., by ATP-driven cross-bridge cycling. To formulate this process as a two-flux-system, a molar torque  m (generated by cross-bridge cycling) is defined as output force A  coupled to the input force A ATPc (A5). An increasing pressure development of the activated heart is considered in simulations through a [Ca 2+ ] c dependent change of  m . From NET it is clear that increasing A  must reduce the resultant force A coup , so that the associated flux must be reduced too. This does not only decrease velocity of contraction, but also increases efficiency of this reaction. Figure 6A shows consumption, which underlines above mentioned results of MCA. With permeabilized cardiomyocytes it could be shown [28] that the activation effect may be caused by a limited availability of ADP for OP, since ADP 3cannot permeate with sufficient rapidity the voltage dependent anion channel (VDAC) [45]. The demonstrated Pi activation of O 2 consumption [46,47] presumably cannot be explained by a limited availability of H 2 PO 4 -, since, in contrast to ADP 3-, this Pi species may easily permeate the VDAC pore of the outer membrane. It was suggested that Pi may activate dehydrogenases of the matrix and in addition may be needed to activate electron transport of the respiratory chain [46]. But even if the above mentioned activation effects of [ADP] c and/or [Pi] c were operative also in coupled glucose oxidation of real VMs, an activation of whole pathway flux is absolutely necessary to deliver free energy under conditions of an increased consumption of free energy. Since it is not stored sufficiently in glucose metabolism itself including OP, it must be delivered by other stores, that is, mainly from the extracellular space in the case of VMs. So the controlled flow of free energy through the whole metabolic pathway is necessary to satisfy changing free energy demands. As is demonstrated with simulations, the most effective points of control are given by [Ca 2+ ] c dependent load reactions and PFK. The flow of information from load to delivery is brought about by the feedback via AK of [AMP] c on PFK.

Forced Oscillations
Up till now all results were related to steady state conditions, intrinsic oscillations do not occur under these conditions and forced oscillations were excluded. VMs in heart tissue, however, contract rhythmically at various frequencies and variable pressure developments. These oscillating contractions are known to be induced by rhythmic depolarizations of the membrane potential ( c ), which in turn trigger periodic Ca 2+ release from and accumulation into the sarcoplasmic reticulum (SR). respectively. Therefore, it seems justified to conclude that deductions made with respect to steady state are valid also under conditions of forced oscillations.

Conductance Matching for Glucose Oxidation and Power Output
Because OP cannot be formulated as a two-flux-system, Equation (10o) would be inappropriate to show that power output of coupled glucose oxidation is near maximal. The problem of non-reliable overall parameters can be circumvented, however, if transformed conductances were taken into account.
It is well known in electrical engineering that power output of a battery becomes maximal when the internal conductance L i matches the conductance of the outer circuit, L e . A battery can be regarded as a coupled two-flux-system, in which a chemical reaction with affinity A 2 is coupled to an electrical potential difference,  = −U, with affinity A 1 . Coupled glucose oxidation, when transformed to equal flux velocities, can be regarded analogously: the input force is given by ov A  , while the coupled output force is given by  . Both forces and conductances i L  and e L  can be obtained from respective dissipation functions as described above. In analogy to a simple electric circuit consisting of a battery whose poles are connected by an electrical resistance (R e = 1/L e ) with an attached load, i L  is related to the whole pathway of glucose oxidation including coupled ATP c production, and L e to ATP c utilizing reactions, i.e., all cytosolic ATPases. As in an electric circuit, only the electrical resistance (R e ) of an energy converter is involved with power maximation. In energy metabolism it is therefore only that part of W L  , which is associated with ATP splitting.
, and using (10d) and (10f), efficiency in terms of  is given by   Figure 9A). In contrast, when PFK activation by [AMP] c is omitted, power output cannot be increased by increasing e L  ( Figure 9C) and, because conductances are not constant, points do not fulfill Equation (13g). Obviously, the feedback by [AMP] c on PFK can activate whole coupled glucose oxidation to appropriately increase ATP production to meet the demands of an increased ATP consumption.  remains remarkably constant at about 0.79 over a nine fold range of power output. It is close to maximal power output. Consequently, also efficiency must be approximately constant and, as follows from Equation (13g), near a value of 0.5 (0.56) ( Figure 9B, D). Interestingly, also efficiency of the glycolytic span from glucose to pyruvate lies near this value.
The Larger deviations of  from 1.0 can be obtained by uncoupling OP. Figure 9F shows that efficiency obeys exactly Equation (13h) over a very wide range of 's. Both conductances, e L  as well as i L  , are increased by uncoupling OP, but e L  to a much higher degree, so that  increases and thus efficiency decreases. Figure 9E shows that also under conditions of markedly increased 's, power output points are close to their hypothetical respective power curve. Comparison of these uncoupled values with those obtained for coupled conditions shows that larger losses of useful power obviously can be prevented. This power preservation is associated, however, with a marked reduction in efficiency ( Figure 9E,F). These results demonstrate in a striking manner, how serious deteriorations of energy coupling can be compensated by an [AMP] c mediated increase of pathway flux. In addition, for the metabolism of VMs they show that maintenance of power output obviously has priority over conservation of efficiency.

Replenishment of [ATP] c
In VMs there are three well known pathways of [ATP] c replenishment: coupled glucose oxidation, which is described above, and two near equilibrium reactions catalyzed by creatine kinase (CK) and adenylate kinase (AK), respectively.
It is widely accepted that the primary tasks of the phosphocreatine (PCr) system is to buffer [ATP] c and to shuttle it between mitochondria and cytosolic locations of high ATP usage like actomyosin ATPases in the myofilament lattice [48,49]. In addition, at least in present simulations, the PCr system is also interconnected with [AMP] c dependent activation of PFK and thus, with activation of whole coupled glucose oxidation. The reduced [ADP] c production through the action of the PCr system counteracts PFK activation. Therefore, delivery of free energy may depend also on this reaction. Perhaps the most important task of the PCr system may be to function as a shuttle for [ATP 4-] and [ADP 3-] c [26], for these compounds do not easily permeate the outer mitochondrial membrane. In addition, they must overcome diffusional restrictions inside the myofilament lattice, where myosin ATPases are located. Reaching of ATPases by diffusion at other locations, e.g., Na-K ATPases at the sarcolemma, or Ca 2+ ATPases at the sarcoplasmic reticulum, may also be aggravated by diffusional restrictions. Since these ATPases are embedded in the filamentous network of the cytoskeleton, which likewise may be of limited accessibility for molecules like ionic ADP and ATP species. Supposedly such structures are charged at their interface between the cytosol and the network and may produce a Donnan potential, which, because of fixed negative net charges, must be negative inside the filament lattice. Because the shuttle molecules may be also charged, any shuttle mechanism has to take into account electrogenic crossing of the outer mitochondrial membrane on the one hand, and charged interfaces of myofibrils and/or cytoskeletons on the other hand. The following scheme describes the principals of PCr shuttling: starting with an increased ATP splitting in the myofilament lattice ('f') by an elevation of [Ca 2+ ] c leads to (in chemical notation) 2 3 In the filament compartment MgATP 2is recycled, whereas PCr 2and Cr are consumed and produced, respectively. In the intermembrane compartment similar reactions occur, but with opposite directions. From the adenine nucleotide exchange reaction (J AE ) at the inner membrane, ATP 4is delivered to the intermembrane space ('im'), and ADP 3is transported into the matrix. Both are recycled by the PCr reaction in this compartment, 4 3 Because of the limited permeability of the VDAC pore for adenine nucleotides, [ATP 4-] im and [ADP 3-] im cannot be equilibrated rapidly with the cytosol. The associated affinity, A ATPim , may be high, so that here the PCr im reaction is forced into the direction shown in reaction (R1d). In contrast, in the environment of ATPases, A ATPf may be lower, whereby here the PCr f reaction may be forced into the reverse direction ((R1a)). At both locations the PCr system may be near equilibrium, but because of different A ATP 's, at different respective mass action ratios in both compartments. ] drives diffusion of this ion from myofilaments to mitochondria. In the intermembrane space, HPO 4 2takes up one H + produced here from the reversed PCr reaction and finally enters as - 2 4 H PO ion H/Pi symport at the inner membrane. As a result ATP, is transported from mitochondria to ATPases and in its split form back to mitochondria. As a prerequisite, in addition to the two-compartment system of classical cellular energetics, a third compartment for A ATP , the mitochondrial intermembrane space is necessary to guarantee gradient formation for PCr shuttling and a high power output. As a result, however, a certain amount of free energy of ATP production becomes inevitably dissipated by the diffusional flows along these gradients.
During repetitive contraction cycles oscillations of respective flows can be expected (Figure 1). It is questionable, however, that [ATP] c and [PCr] oscillations occurred with such high peak values of about 10% as was found by Honda et al. [50], as this would require extreme rates of O 2 consumption during recovery phases. When the metabolism of VMs is changing from a low to a high power output, the rate of ATP c consumption becomes elevated above the rate of ATP c production, so that A ATPc is lowered, and the PCr reaction begins to partially compensate the decline of A ATPc . After a certain delay, consumption and production rates match again, and A ATPc is adjusted to a new lowered value, which forces the corresponding near equilibrium PCr reaction to a new mass action ratio with a lower [PCr]. A CPK oscillates now again around zero. When power output switches back to lower values, the reverse processes lead to a re-increase of A ATPc until respective ATP c rates match. The PCr reaction again oscillates around zero driving force, but now at a re-increased [PCr]. Supposedly, it is not necessary that concentrations of the PCr system must be changed during a contraction cycle to allow formation of sufficiently steep gradients. As mentioned above, this might be brought about by compartmentalization of A ATPc .
Gradients of [PCr 2-] and [Cr] may be produced preferably at interfaces of mitochondria/cytosol and myofibrils/cytosol, respectively. For the maintenance of flows, it is necessary that compartments of high and low A ATP are in close proximity, because otherwise gradients could be destroyed by soluble CK of the cytosol which may be present also in this compartment at high activity. Such structural prerequisites obviously are realized in VMs, in which mitochondria surround myofibrils and the sarcolemma in a highly ordered way.

Pancreatic β-Cells
The primary task of β-cells is to secrete insulin in response to an elevated glucose concentration. The process of secretion involves exocytosis of insulin containing vesicles, which are [Ca 2+ ] c and [ATP] c dependent like other load reactions of the cytosol. The sequence of events, however, leading to an activating increase of [Ca 2+ ] c in β-cells is remarkably different from activation of VMs. In contrast to this latter cell type, an increase of [Ca 2+ ] c in the β-cell is induced intracellularly by an increase of pathway fluxes of coupled glucose oxidation and associated decrease of [ADP] c (increase of A ATPc ). The lowered [ADP] c in turn closes ATP-sensitive K + channels (K ATP ), leading to a depolarization of  c  and an influx of Ca 2+ from the extracellular space through opened (by depolarization) voltagegated Ca 2+ channels (Ca V ) into the cytosol [51]. Such a metabolically induced mechanism of [Ca 2+ ] c increase can certainly only be realized by oscillations. A low [ADP] c and an elevated [Ca 2+ ] c cannot coexist over a longer period, because an increased [Ca 2+ ] c would activate cytosolic load reactions including insulin exocytosis, which in turn leads to a re-increase of [ADP] c , a repolarization of  c  via opening of K ATP channels, and a re-closure of Ca V channels. A metastable state must be generated, which is characterized by oscillations with a phase shift between the two variables [Ca 2+ ] c and [ADP] c [6,10]. That is to say, the question why stimulus-secretion coupling of the β-cell must be associated with oscillations is explained by the metabolically induced reaction sequence, which itself represents the intrinsic part of an oscillator.
The initiation of an increased glucose metabolism, in contrast to VMs, does not proceed via activation of PFK by [AMP] c , but through activation of glucokinase by cytosolic [Glu]. Glucokinase has a much higher K M value than hexokinase from heart muscle and, therefore, is not saturated at resting acts as a stimulus for insulin secretion, whereas the glucokinase reaction, beside its function as a hexokinase in GLY, constitutes in addition the recognition mechanism for this special stimulus.
Glucokinase of β-cells is known to be not inhibited like hexokinase by [G6P]. In simulations of glucose metabolism of VMs, the inhibition by the latter variable is needed to make this first reaction step more flexible for the adjustment of a steady state. Because [Glu] is not a variable but a fixed parameter, flexibility of the first step of GLY in β-cells has been introduced by contracting the first three reactions of GLY, including the PFK reaction to one single step (see [6]). Now the first step of GLY in simulations possesses an activation factor containing [ADP] c as a variable, so that certain flexibility also for β-cell simulations may be ensured. This is the main difference (beside different geometrical parameters and conductances) between simulations of whole coupled glucose oxidation of VMs and β-cells. In these latter cells the concentration of total adenine nucleotide concentrations is, however, much lower than in VMs. In β-cells, especially [ATP] c is much lower (in simulations 3,740 µM compared to 8,820 µM). The responsiveness of PFK and PK to the inhibitory action of [ATP] c , therefore, may be drastically lowered in β-cells, even if in this cell type the same isoenzymes as in VMs were expressed. That is, the functional characteristics of these enzymes might be markedly changed simply by a lowering of [ATP] c .
The above described mode of glucose recognition has its price, however. Compared with VMs, delivery of free energy in β-cells via coupled glucose oxidation is drastically cut back through the limited capability of GK to respond on During periods of hypoglycaemia, β-cells would be even more easily affected. To counteract such a risk, it is indispensable that delivery of free energy can be partly provided by a parallel fuelling of other substrates like fatty acids and/or amino acids into energy metabolism.
At a [Glu] of 10.0 mM β-cells are less susceptible to cell damage, but insulin release can be seriously jeopardized by a relative weak uncoupling. Under activated conditions at very low values of  sp =  sh = 0.015, oscillations disappear and insulin exocytosis is decreased to 7.6% (not shown).
From this behavior of β-cells, a fundamental conclusion with respect to the degree of coupling of mitochondrial inner membrane reactions can be reached: especially ATP synthase must be tightly coupled, because otherwise [ADP] c in β-cells would be markedly increased. Furthermore, because it seems plausible to suggest that mitochondria of other cell types like muscle cells might not differ in this respect, a tight coupling of inner mitochondrial membranes can be expected.
As is already described above, the PCr system shuttles ATP free energy between mitochondria and ATPases. It is questionable, however, if this system of facilitated transport may be operative in β-cells, because the buffering effect of the PCr system might interfere with the oscillatory machinery of stimulus-secretion coupling. Figure 10 shows results from a β-cell simulation (taken from [10]), in which the PCr system additionally has been included. Stimulation of β-cell activity by 10.0 mM [Glu] is possible only if [PCr tot ] does not exceed 1.0 mM. At higher concentrations normal oscillations of  c  and [Ca 2+ ] c cannot be elicited anymore, whereby insulin release becomes markedly reduced. In the absence of PCr shuttling, normal functioning of stimulus-secretion coupling nevertheless may be possible, since the oscillatory frequency of activated β-cells is very low compared to rhythmic contractions of VMs, so that diffusional restrictions for adenine nucleotides cannot impair power output.

Methods
For a treatment of coupled glucose oxidation as an integral part of cellular energetics, cell membrane and SR reactions are of secondary importance. Therefore, to simplify expressions, all reactions of the cell membrane and  c  are omitted in present simulations. In addition, also most of SR reactions are excluded, only Ca 2+ pumping of the sarco/endoplasmatic reticulum Ca 2+ ATPase (SERCA) as an ATP consuming reaction is included in simulations. [Ca 2+ ] c therefore appears in calculations not as a variable but as an adjustable parameter.
Because of different geometrical parameters, volume fractions and  values have to be adapted for VMs. Data from Aliev et al. [52] are used to calculate volumes of cellular compartments. A value of 50 and 20% of V Cell (V Cell = 34.4 pL, [53]) is used for sarcosolic volume V c , and for total mitochondrial matrix volume V m , respectively. Capacitances of all mitochondrial inner membranes of 3.6 nF and   Fluxes J AK and J CPK (A10 and A11) are added to a simulation taken from reference [10].

Conclusions
Applied equations are structurally similar to other phenomenological relations and especially to equations used in NET. But in contrast to these latter equations, the new equations contain variable conductances, which allow modeling of individual fluxes in analogy to the kinetic behavior known from enzyme-catalyzed reactions. At steady state, all variables are constant, so that flux equations presented here become undistinguishable from those of NET and Michaelis-Menten type equations. The great advantage over the latter rate laws, however, is given by its remarkable stability towards perturbations and its ability to reach a new steady state after such a perturbation in a possible short interval. This behavior of equations makes them suitable for modeling studies of metabolic networks and, in addition, for elucidation of pathway control mechanisms.
In simulations of coupled glucose oxidation, a drastic increase of free energy flow to increase power output can be achieved only at two points, at [Ca 2+ ] c dependent load reactions and at PFK. For maximal power output respective conductances of in-and output fluxes must be similar in magnitude. This conductance matching is sufficiently fulfilled over a wide range of useful power generation. It may represent the basic mechanism of energy metabolism to ensure a high cellular power output under conditions of a sufficient supply of substrate and oxygen. Moreover, it may be regarded as a necessary prerequisite for the development of so-called fight or flight behavior, as it can be observed in higher vertebrates as well in more primitive forms.
In β-cells, the GK enzyme is limiting the metabolic rate of coupled glucose oxidation, whereby the metabolic mechanism of glucose recognition is made possible. Obviously, only a few changes in the first reaction steps of glucose metabolism are sufficient to produce that remarkably large difference of power output between VMs and β-cells. On the other hand, both cell types differ enormously not only in size and geometry, but also in their markedly different degrees of structural organization. This strongly indicates that a task such as high power output may be guaranteed only if, in addition to adequately controlled pathways, that appropriate structural realities also exist. Only in this way, conductance matching for a maximal power output can be realized.

Ventricular Muscle Cells
Nearly all individual fluxes of coupled glucose oxidation are given in references [6] and [10]. In VMs only the PFK reaction is added to these pathway fluxes, leading also to changes of the formerly more contracted flux through GK plus PFK.     The above flux equations show that often several activating or inhibiting factors may be necessary to model a given flux. Principally these probability factors can be replaced by one single factor ( max Enz v V  , v = J M-M ) obtained from the given rate law of enzyme kinetics.

Oxidative Phosphorylation
The fraction Q H of input fluxes J NA plus J FA for OP is given by

OP as a Two-Flux-System
Conductances L OP and L Wop and overall  values (1 ov and (2 ov , respectively) can be obtained from dissipation functions. OP then can be formulated as a two-